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TECHNICAL MEMORANDUM 


MODELING OF TURBULENCE EFFECT ON LIQUID JET ATOMIZATION 

1. INTRODUCTION 


Recent experimental investigations and physical modeling studies have indicated that turbulence 
behaviors within a liquid jet have considerable effects on the atomization process. Such turbulent flow 
phenomena are encountered in most practical applications of common liquid spray devices. Most of the 
existing atomization models do not account for turbulence effects. Limited attempts have been made 
to model the subject effects on liquid jet disintegration; however, they treat turbulence as either an only 
source or as a primary driver in the breakup process. 

This research aims to model the effects of turbulence occurring inside a cylindrical liquid jet 
to its atomization process. Although considerable research has been oriented to liquid atomization 
processes for more than a century by many authors, only work that is closely relevant to this investiga- 
tion is considered and presented here. In summary, this research effort contributes to the improvement 
of predictions in liquid jet atomization processes. 

Outlines of the proposed research topic will be presented in the following subsections. The 
motivation and objective of this research effort are stated along with a discussion of the approach to 
be taken. 


1.1 Background 

The transformation of a liquid body into droplet sprays in a gaseous surrounding is of great 
importance in many industrial processes. Various spray devices have been developed that are designated 
as atomizers or injectors. In liquid fuel combustion systems, such as rocket engines, diesel engines, gas 
turbines, and industrial furnaces, the combustion efficiency and chemical reaction behavior are primarily 
dependent on the effectiveness of the liquid body broken up into sprays. Smaller drop sizes, generated 
from the spray devices, increase the specific surface area of the fuel and thereby achieve high rates of 
evaporation and mixing. Consequently, higher engine efficiency would be realized in the combustion 
process. On the other hand, mixed fuel that reacts rapidly near the injector could cause the injector exit 
surface to inadvertently overheat. These phenomena have been observed in many liquid rocket engines 
during their hardware development phase. Therefore, understanding and adequately predicting this 
physical breakup process leads to designing better spray devices for these various applications. 

In general, the transformation of a liquid jet into droplet sprays can be divided into two separate 
stages as portrayed in figure 1. Lefebvre describes the mechanisms of the jet breakup as follows: “When 
a liquid jet emerges from a nozzle as a continuous body of cylindrical form, the competition set up on 



Primary Liquid Jet Breakup 



Secondary Drop 
Breakup 


Figure 1. Sketch of primary and secondary liquid jet breakup regimes. 


the surface of the jet between the cohesive and disruptive forces gives rise to oscillations and perturba- 
tions. Under favorable conditions the oscillations are amplified and the liquid body disintegrates into 
drops.” 1 This process is sometimes referred to as primary atomization or primary liquid jet breakup. 
Moving downstream, the liquid further breaks up into smaller drops. This process is also known as 
secondary atomization (secondary droplet breakup). 

The mechanisms of atomizing the liquid jet are complex. Jet inertial and aerodynamic forces 
along with the surface tension certainly play a role in the breakup process. In addition, interaction of the 
liquid with turbulent gas flow also contributes to the jet disintegration. Considerable studies have been 
dedicated to gain knowledge and subsequently model this phenomenon. Faeth and Crowe et al. provide 
substantial reviews on this topic. 2 ' 3 

Many researchers have also reported that for certain favorable flow conditions turbulent motion 
can be recognized within the liquid jet. Often the geometrical sharpness of the injection nozzle inlet, 
along with appropriate flow conditions, can create cavitations inside the nozzle. The collapse of this 
cavitation can generate a flow fluctuation, leading to a more aggressive disintegration of the liquid jet. 
The turbulent primary breakup was first reported in 1932 based on visualization studies conducted by 
De Juhasz et al. 4 In the following years, Schweitzer observed that the level of the liquid turbulence had 
some effect on the jet stability. 5 Later investigations of the turbulent liquid jet breakup in quiescent air 
were performed by various authors, including Chen and Davis, Grant and Middleman, Phinney, and 
McCarthy and Mallory. 6-9 They observed that disintegration of the liquid jet occurred rapidly near the 
nozzle exit in certain flow regimes. Using the measured data, Grant and Middleman were able to corre- 
late the flow conditions and liquid breakup lengths. Phinney conducted a series of experiments in which 
the liquid jets were propelled into a gaseous atmosphere from injection tubes at various fluid conditions 
ranging from laminar to fully developed turbulent flows. Similar to Grant and Middleman, Phinney also 
came up with a correlation of breakup lengths based on dimensional analyses. His results showed that 
the breakup length of the liquid core was significantly shortened when the liquid jet changed from 
laminar to turbulent flow. The author also observed the appearance of fine milk-white droplet clouds at 
the injection ports. 

With the advancements of nonintrusive measurement techniques, introduced in the 1980s, 
researchers have been able to experimentally discover more physical details and provide more measured 
data for subject atomization. Studies of relatively large-size liquid jets, performed by Wu et al. and 
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Tseng et al., provided droplet data, such as size, velocity, and fluctuation quantities, and reported a more 
complete mechanistic approach for the primary atomization regime. 10-14 During the investigation of fully 
developed turbulent round jet sprays, Wu et al. and Tseng et al., observed that when turbulence 
developed inside the jet column, starting at the nozzle exit, this turbulence motion remained dominant 
and became a main contributor in the spray development. The authors believed that this turbulence 
characteristic played a primary role on the liquid stripping near the injector exit. Tseng et al. also 
indicated that drop sizes after the primary breakup, in general, were larger for turbulent than nonturbu- 
lent liquids. Based on the data collected from their experiments, the authors were able to establish 
correlations between the turbulence fluctuation quantities, the breakup drop size, and the breakup length 
of the liquid jet. Further studies of turbulent liquid jets flowing into a quiescent gas were also performed 
by Sallam et al. 15,16 They characterized the liquid breakup by a series of flow visualizations as well as 
phenomenological analysis. The results showed that the breakup regimes for the turbulent liquid jets can 
be described in a way similar to the existing theories of nonturbulent liquid jets, and these regimes can 
be fully expressed by the turbulence flow properties. 

1.2 Previous Research on Atomization Modeling 

In the arena of atomization modeling, various theories and models for cylindrical jet atomization 
have been proposed. In 1879, Rayleigh studied the linear stability of an infinitely long, circular, inviscid, 
incompressible liquid jet. 17 He was able to derive an expression for the growth of surface disturbances as 
a function of time. This theory works reasonably well for the low-Weber number regime. Theoretical 
work and model development, based on the Kelvin-Helmholtz (KH) instability of the liquid surface 
wave, has continued to grow significantly over the years. Comprehensive reviews on this classical work 
can be found in McCarthy and Bogy. 9 ' 18 Later, Rangel and Siringnano, and Clark and Dombrowski 
extended the surface instability theory for the liquid jet and sheet breakup to include higher order terms 
that account for the effects of surface tension and liquid-to-gas density ratio. 19-21 

For the last few decades, developments and enhancements of ink-jet printing, diesel fuel injec- 
tion, and direct injection techniques in automobile engines and liquid rocket engines have required 
further improvements of liquid injection predictions. In addition, improved computational techniques 
and advanced computer capabilities have promoted further research in this modeling arena. As to the 
convenience of implementing improvements to computational fluid dynamic (CFD) methods and the 
wide utilization by analysts, the two most noticeable atomization models are the KH instability of Reitz 
and the Taylor analogy breakup (TAB) model of O’Rourke and Amsden. 22 2 ’ Reitz derived the KH 
model, also known as the blob model or stripping rate model, which described the primary breakup 
entirely in terms of the surface wave perturbation, in which viscosity, surface tension, and aerodynamic 
forces were the contributing factors. On the other hand, the TAB model was based on an analogy 
between an oscillating distorting droplet and a spring-mass system. In this model the drop distortion 
was driven by the force interaction between the external aerodynamic surface tension and the viscous 
damping of the droplet liquid. The TAB model is suitable for predicting the secondary breakup regime. 
Further reviews of these two models can be found in section 2. Recently, Tanner examined data from 
liquid jet breakup experiments at high-pressure conditions and found that an intact liquid core is broken 
into various liquid shapes or drop sizes shortly after the injection exit port. 24 Hence, the author extended 
the TAB model and developed an enhanced TAB (ETAB) model to include the primary breakup regime. 
Tanner also accounted for the droplet surface stripping near the injector exit with a power law drop-size 
distribution using his recent cascade atomization breakup (CAB) model for high-pressure liquid jets. 25 
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Though these models in some aspects provide reasonable predictions of the liquid atomization, 
they do not account for the liquid turbulence motion observed in certain sprays as previously described. 
In a study of primary diesel fuel atomization, Nishimura and Assanis developed a phenomenological 
cavitation model in which the primary breakup is governed by the turbulence mechanism . 26 Turbulence 
energy is formed by bubble collapse and fluid turbulence motion. To consider the nozzle exit turbulence 
conditions of diesel spray modeling in the atomization process, Huh et al. proposed a scheme taking into 
account two independent mechanisms— wave growth and turbulence . 27 The turbulence is characterized 
partially by the injection nozzle geometry, while the wave growth is derived from the KH instability 
theory. A logical framework of coupling the flow inside the injection nozzle to the jet breakup can be 
achieved with this approach. The rationale for this approach is based on an order of magnitude analysis 
of the flow dynamic breakup mechanisms and those associated with the turbulence. This analysis 
concludes that gas inertia force and turbulent stress are the two main forces in the atomization of typical 
diesel engine injectors. In this model, the breakup rate for the parent drops is set to be proportional to the 
length scale-to-time scale ratio. It is hypothesized that the turbulence length scale is more dominant in 
the primary breakup while the wavelength scale is more important for the secondary droplet breakup. 

On the other hand, the time scale is a linear sum of the turbulence and wave growth time scales. 

The aforementioned models were formulated with semiempiricism. The associated coefficients 
had to be determined with reference to experimental data. Another class of atomization models was 
derived based on the conservation laws of mass, momentum, and energy. By taking advantage of 
improvements in CFD techniques, as well as new computer capabilities, the hydrodynamic equations 
that account for nonlinear effects can be directly solved. Various numerical computations using the 
volume-of-fluid method, as proposed by Fromm, Child and Mansour, Shokoohi and Elrol, and more 
recently Liang and Ungewitter, represent the first attempts of direct simulation without using the wave 
instability theories . 28 31 In their analyses, the Navier-Stokes equation system has been numerically solved 
to predict the interfacial instability and jet breakup. Their results generally agree with the measured data 
for low-velocity injection conditions. Yang et al. proposed to solve the jet hydrodynamic equations 
which also accounted for the nonlinear inertia terms and full nonlinear capillary pressure terms . 32 

Several efforts using more sophisticated numerical approaches have been performed recently 
modeling the detailed turbulent flow fields in liquid and gas during the atomization process. Klein 
And Janicka and De Villiers et al. computed the two-phase flows using the large eddy simulation (LES) 
method . 33 34 For the first time, De Villiers et al. were able to apply a numerical technique to resolve jet 
atomization under the more realistic operating conditions found in diesel engines. Leboissetier and 
Zaleoki also attempted to simulate a liquid spray with a multidimensional pseudo direct numerical 
simulation (DNS) method . 35 

The aforementioned numerical methods have a real potential of providing a complete physical 
description of the liquid jet breakup with minimum assumptions; however, they require small computa- 
tional time steps and fine grids across the entire jet domain for their simulations. The LES and DNS 
techniques, in particular, may need submicron spatial elements in size and pico-second in time steps 
to properly predict the atomizing sprays at the high-velocity injection conditions. Consequently, grid 
mesh size and the considered physical domain must be taken care of so that the computational time and 
memory storage requirements will be manageable. At the present time these techniques are still too 
expensive and generally impractical in terms of computational time and power requirements for 
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engineering calculation applications. Hence, engineering analysis and design of the liquid spray devices 
must still rely on physical engineering models. 

1.3 Motivation and Objective 

The previous reviews demonstrate that researchers have observed the natural existence of turbu- 
lence within liquid jets for certain flow conditions. Generally, this flow fluctuation can be formed by 
favorable flow conditions, for instance, a fully developed turbulent flow in a long cylindrical tube. 
Sometimes it can also be created from the collapse of cavitation bubbles inside the injection nozzle. The 
cavitation is induced by the shape of the nozzle, such as a sharp nozzle entrance or nozzle curvature. 
Unfortunately, the aforementioned turbulent flow conditions are present in most applications of common 
liquid spray devices. Moreover, results of past experiments suggest that turbulence motion within the 
liquid may alter the intact jet core, as well as the production of drops that have a different size in com- 
parison to the droplet breakup without the turbulence consideration. Undoubtedly, an accurate modeling 
of such turbulence effects on the atomization process would significantly enhance the prediction capa- 
bility of existing physical models. Limited attempts have been made in modeling the turbulence phe- 
nomenon on liquid jet disintegration. Nonetheless, subject correlation and models treat the turbulence 
either as an only source or as a primary driver in the breakup process. Other computations of the liquid 
atomization have been performed using DNS; however, these sophisticated numerical simulations are 
computationally unaffordable at the present time. 

Hence, the present research aims to model the turbulence effects on the primary and secondary 
atomization processes. Such models can be implemented in numerical simulations to predict the breakup 
process of a cylindrical liquid jet. This proposed effort enhances the existing atomization models so they 
can predict the liquid breakup with more physically realistic conditions. 

1.4 Research Approach 

In the course of this study, terms accounting for the turbulence motion within a liquid were 
developed based on a phenomenological point of view. These terms are appropriately supplemented 
to the two classical atomization models— KH primary breakup of Reitz and TAB secondary breakup of 
O’Rourke and Amsden . 22 ' 23 In the primary atomization model, the level of turbulence effects, on the 
mass stripping of blob drops and product drop sizes are represented by the characteristic time and length 
scales and kinetic energy. This treatment offers contributions of individual physical phenomena on the 
liquid breakup. For the secondary breakup, an additional turbulence effect, acted on the parent drops, is 
modeled and integrated into the TAB governing equation. This turbulence term is referenced to the 
dissipation rate of the turbulent kinetic energy and the deformation rate of the parent drop distortion 
displacement. All associated empirical constants used in these additional terms are determined by 
matching the results with available measurements. Two computer programs have been written to simu- 
late these proposed models. Results of the new and existing models are compared and fully appraised. 
Finally, these models are implemented into existing CFD simulations. The predicted results are evalu- 
ated with the available experimental data. 
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2. REVIEW OF EXISTING ATOMIZATION MODELS OF LIQUID JETS 


This section provides an extensive review of the KH instability model of Reitz and the TAB 
model of O’Rourke and Amsden, since the proposed models are based on these two methods . 22 21 The 
material in this section provides complete descriptions and formulations of the two subject models. The 
governing equations of the liquid surface wave motion are reviewed leading to the blob primary jet 
breakup model. Expressions for the mass stripped from the blob drop and the produced drop size are 
also presented. The TAB model derivation is reviewed, concentrating on the forces acting on a liquid 
droplet immersed in a gaseous flow field. The equation for estimating the product drop size is also 
formulated. As stated earlier, the KH instability and TAB for primary and secondary atomization are the 
basic foundation for formulating the turbulence terms included in these models. 

2.1 Kelvin- Helmholtz Instability: Primary Atomization Model 

The original formulation of this primary breakup model was based on the stability analysis per- 
formed by Reitz and Bracco . 36 They examined the surface stability of a cylindrical liquid jet that is 
subject to perturbations using a first-order linear theory. The surface perturbation and consequential drop 
breakup can be visualized by the sketch shown in figure 2. This sketch shows a blob of liquid with 
surface waves at its interface. The waves continue to grow in amplitude until the breakup occurs and 
then the product drops are formed. The analysis starts by imposing an infinitesimal axisymmetric 
perturbation displacement of the form on the surface: 


r 1 = 9t(r ]0 e ikx+a,t] j , 


( 1 ) 


where 


i=yf— I . 


The displacement ( rj ) is a real part (,9i ) of a function, which contains variables of the time (/) and the 
axial distance (x) in reference to the initial jet injection time and respective nozzle exit location, rj is 
related to the growth rate (a>) of an initial perturbation of infinitesimal amplitude ( 770 ) to its wave 
number ( k ). 


By applying the linearized hydrodynamical equations for the liquid and gas, the inviscid gas 
equations of motion yield for the gas pressure at the interface r=a 


(O 


Pg = ~Pg W -i — kr] 


Kq ( ka ) 
K[ (ka) 


( 2 ) 
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Figure 2. Schematic showing surface waves and breakup of a liquid blob. 


where r is the radial coordinate of the cylindrical liquid jet, while p g and are the gas pressure and 
density, respectively. W is the magnitude of the relative velocity between the liquid and gas. K {) and K\ 
are zero and first-order modified Bessel functions of the second kind. With the assumption that ?/ « a, 
the kinematic, tangential stress, and normal stress equations at the interface can be expressed as 

d ?7 dm dvj 

Vi = — , = (. 

dt dr dx 


and 


~Pl + 2v lPl 


dvj_ 

dr 


O 

7 


?7 + a 


2 d~?7^ 
dx 2 


Pg 


= 0 . 


( 4 ) 


The axial and radial liquid velocity components u/ and V[ are along the liquid surface. Surface 
tension, dynamic viscosity, density, and pressure of the liquid are a, 1 7 , pj, and pj, respectively. 

Solving equation (4) along with the substitutions of equations ( 1 )— (3) results in a dispersion relationship 
as follows: 


7 


co 2 + 2v,k 2 co 


7j {ka) 2 k$ I\ {ka) I[ (ksa) 


Iq {ka) k 2 + ks 2 Iq {ka) 7, {3a) 


ok 

Pl a *~ 

p S ( , 


i 1 2 2 

1 - k a 


!s 2 - k 2 ) 

I\ {ka) 

^ks 2 + k 2 j 

h) {ka) 

,2^ 2 - 

k 2 ) I Y (ka) K { {ka) 

1 3 2 + k 2 Io{ka) K 0 {ka ) 


( 5 ) 


where 


3 = \Jk~ +co/v, . 

The above equation contains 7 0 and 7 1 as the zero and first-order modified Bessel functions of the first 
kind. The first derivative with respect to time for 7 1 is 7{ . Reitz numerically solved equation (5). 22 The 
results indicate that there is a single maximum in the wave growth rate curve co= co (k). He further 
generated curve fits of the numerical solutions to equation (5) for the maximum growth rate (co = Q) 
and for the corresponding wavelength (A = A) that led to the correlation as follows: 


— = 9.02 
a 


1 + 0.45Z 05 111 + 0.4T 


, 0.7 


1 + 0.87We 


1.67 
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and 
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Pia 3 

03 | c 
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where 


Z 

Re } 

We, 

We g 

T 


Ohnesorge number \ Z = We/ 0 ' 5 / Re, 

Reynolds number of the liquid {Re, = Wa/v ,) 

Weber number using the liquid density [we, = p,W 2 a/o 
Weber number using the gas density [we g = p g W 2 ajo 
Taylor parameter i'ZWe^ 


(6) 


(7) 
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To implement equations (6) and (7) in the primary atomization prediction, Reitz assumed that 
the liquid was injected in a monodispersed fashion, of which the liquid jet initially was in a form of blob 
parcels, containing liquid spherical drops with their size equal to the injection nozzle diameter. Hence, 
the Reitz model was also interchangeably called the blob model. 

Traveling downstream from the injector, the blob drop (also referred to as the parent drop) 
undergoes surface breakup into product drops (also denoted as children drops) due to a tripping action. 
This is due to the unstable surface wave driven primarily by the fastest wave growth rate and the corre- 
sponding wavelength. The rate of change in the parent drop size is governed by a ratio of the character- 
istic length and time scales associated with the wave perturbation motion as 

da a - r,„ 

— = , when r w <sc a . (8) 

dt r 


The time (r) is associated with the surface wave motion while r w is the radius of the product droplets 
that are formed from the shredded mass of the parent drop. The parameters r w and r are related to the 
fastest wave growth rate and the corresponding wavelength as 


[B 0 A 


min 


3jra 2 w/2n)° 33 , 


( BqA < a) 
0.33 


( BqA > a) 


(9) 


and 


x = 3.126B x a/AQ , (10) 

where the values of Bq and B\ are set equal to 0.61 and 10, respectively. The approach used in this study 
is directed to extend this formulation by including a turbulence term, as shown in the next section. For 
that reason, it is helpful to rewrite equation (8) as follows: 


da 

dt 




, when BqA a . 


( 11 ) 


L w = A and r w = aj AQ can be considered as the characteristic length and time scales associated with 
the unstable surface wave behavior. The coefficient C a becomes 


a 3 . 726 #! 

In summary, the blob model that is used to describe the primary breakup was formulated entirely 
in terms of the surface wave perturbation in which viscosity, surface tension, and aerodynamic forces 
were the contributing factors. Hence, the present form of this model needs to be modified when turbu- 
lence appears in the liquid jet. The proposed enhancement will be discussed in section 3. 
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2.2 Taylor Analogy Breakup: Secondary Atomization Model 

O’Rourke and Amsden constructed the TAB model based on an analogy between an oscillating, 
distorting droplet and a spring-mass system. Formulation of the subject models is presented in the 
remainder of this section. 

2.2.1 Modeling Droplet Surface Wave Motion 

In the TAB model, due to the surface tension and viscous damping of the liquid droplet, drop 
distortion is driven by external aerodynamic droplet gas interaction force and droplet shape restoring 
force. The analogy between the droplet distortion and a spring-mass system is shown in figure 3. 




Figure 3. Schematic showing analogy of forces acting on a liquid drop and spring-mass system. 


The linear differential equation for a forced dampened harmonic oscillator is written as 

= F - k£- dt , (12) 


where 

m = mass of drop 

£ = radial cross-section change from its equilibrium position (displacement of the equator 
of the droplet from its equilibrium position) 

F = force acting on the drop due to the aerodynamic droplet-gas interaction 
k'C = droplet shape restoring force caused by the surface tension 

ciC, = damping force created by the liquid viscous effect. 
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Formulations for the coefficients used in the terms of equation (12) are 


F 

m 



Pl r p 


k 

m 


C k 


a 

Pl r p 


d r n 

~ = L d 9 

m Pl r p 


where and pj are the gas and liquid densities, W is the relative velocity between the local gas and 
the droplet, r p is the droplet radius, a is the liquid surface tension coefficient, and Uj is the liquid 
viscosity. Based on matching the analytical results with experimental data the constants C F , Q., and Cp 
are set to 1/3, 8 , and 5, respectively. 

To nondimensionalize the distortion displacement (£), y is set equal to , where C/, is 

a constant. When substituting this nondimensional parameter and the coefficient formulations into 
equation ( 12 ), it can be rearranged as 


>’ = 


C F P g W z 
c b Pi rl 
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~ c k — 3 y~ c d 

Pi r p 


pi 

Pi r p 


■y ■ 


(13) 


Equation (13) is a linear, nonhomogeneous second-order differential equation. It has an exact solution as 

>’0 
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(14) 
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where 


PnW 2 r p 

We = ^ l - 

o 

1 _ C d Hi 
l d 2 Pl>'p 

f 2_ r ° 1 
0) Cfc . 

pi>'p t d 

The Weber number and angular frequency of the oscillatory motion are denoted as We andw . The 
values of y and its first derivative at / = 0 are vq and y 0 . 

2.2.2 Drop-Size Estimation From TAB Model 

The Sauter-mean radius (SMR) of the product drops is derived from the energy balance before 
and after droplet breakup. Prior to the breakup, the energy of the parent drop is represented by the 
minimum surface tension energy (E sur j) and the oscillation and distortion energy ( E nsc ) that are defined 
as follows: 


E surf 4 ^ r p & ( 15 ) 

E„ sc - +m 2 C 2) j- K^pif [y 2 +<o 2 y 2) j . (16) 

The constant (K) is the ratio of the total energy from distortion and oscillation to the energy in the 
fundamental model. The value of £ = 10/3 has been obtained by matching analysis data with experiment 
data. So, the total energy of the parent drop prior to the breakup (E par ) is 

Epar = E swf + E osc = 4 rfo + K^p,r p (y 2 + co 2 y 2 ) . (17) 

It is assumed that the product drops are not distorted or oscillating after breakup. Thus, the 
energy after breakup (E prod ) would only contain the minimum surface energy and the kinetic energy 
of the motion normal to their parent droplet 


E prod 4jzr p O 


r 32 


TC 5 .7 

+ ~ r P piy 


(18) 
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By equating equations (17) and (18) the SMR (r 32 ) of the product drops, formed after the secondary 
atomization, can now be calculated as 


77 _ JO 

^ prod ^ par 

4:tr„o 

r 32 = . (19) 

rr 71 5 .2 

E par g r p p0 

As mentioned in section 1, the liquid turbulence phenomenon has been encountered in certain 
spray conditions. Earlier investigations of such turbulence characteristics suggest that these turbulence 
behaviors also contribute to the jet breakup process. 81013-16 Fluctuation of the liquid properties can be 
generated from the turbulent flow conditions that exist inside the injection nozzle. In addition, injection 
nozzle hardware geometry and surface roughness can increase the fluctuation; however, the present 
analytical formulations of these models are not able to capture this physical behavior. This research is 
focused to model the liquid turbulence effects and incorporate them into existing breakup models, with 
the goal of improving the breakup process prediction in a more realistic fashion. Modifications of the 
KH and TAB models are proposed in section 3. 
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3. MODELING OF TURBULENCE EFFECT ON LIQUID JET BREAKUP 


This section presents the derivation of the primary atomization model (Tblob) for the turbulent 
liquid jet. All physical models and closure expressions associated with this breakup regime are shown. 
This section focuses on formulation of the terms to account for the turbulence effect on the primary 
atomization. These terms are derived based on the phenomenological approach. The expression repre- 
senting the rate of change in the blob drop size is modified and the turbulence characteristic scales to be 
used in the primary breakup model are formulated. The proposed estimation of the product drop size and 
its velocity are also presented. Finally, the expression for estimating the initial turbulence quantities of 
the product drop is developed. 


3.1 Rate of Change in Parent Drop Size 

The main contribution of this research effort is oriented to incorporate the turbulence effect in 
modeling of the liquid jet breakup. As mentioned in section 2, the blob model is extended to include 
terms and parameters associated with the turbulence behavior. Hence, it may be proper to designate the 
present model as the Tblob model. The phenomenological approach is chosen in this model develop- 
ment to account for both the effects of surface wave perturbation and turbulence motion. The resulting 
model should adequately include the combination of these two phenomena when both coexist in the 
considered flow conditions; however, the model is also capable of describing the breakup caused by 
surface wave perturbation when the turbulence is not present in certain flow conditions. In such a case, 
the original formulation should be retained. For all aspects, the model should reasonably maintain a 
smooth transition from the nonturbulent to fully turbulent flow regimes. Also in this phenomenological 
method, the kinetic energies associated with surface wave and turbulence motions are used to weight 
these individual effects. Consideration of the kinetic energy level for this weighting is based on the 
phenomenon that the larger kinetic energy motion would have a stronger influence in the liquid jet 
breakup process. 

In the proposed primary breakup model, length and time scales are composed of those described 
in the KH model, representing surface wave instability and turbulence behavior. Subject turbulence is 
characterized using the framework developed by Huh et al. 27 

To derive the new model, it is proposed that the rate of change in the parent drop radius has a 
similar formulation as shown in the Reitz model, equation (11), with an inclusion of the turbulence term 


da 

dt 


--C 

( T 




T 


T t ). 


, when r p «; a 


( 20 ) 


L t and r t are the turbulence characteristic length and time scales, respectively, and r p is the radius of the 
product drop. These parameters will be defined later in this section. The third term on the right-hand 
side of equation (20) represents the contribution of the turbulence effect on the stripping rate of the 
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parent drop. This turbulence term would accelerate the reduction in the parent drop radius with time. 

As expected, the primary breakup should be complete earlier than for the case without the turbulence. 
This result is consistent with the experimental observations of Phinney, Soteriou et al., Chaves et al., and 
Schimdt et al. in which the liquid jet has a short intact core length when turbulence appears due to the 
flow at a high Reynolds number or cavitation. 8,37-39 The contribution of the unstable surface wave to the 
primary breakup is reflected in the first and second terms on the right-hand side of the equation. They 
are the same terms as shown in the Reitz model. It should be noted that the third term would vanish 
when turbulence does not exist in the fluid. In such a case, equation (20) would become the original blob 
model of Reitz that only retains the effect of the surface wave motion on the primary breakup. 

Similar to the criterion set by Reitz, the parent drop would no longer strip its mass to create the 
product drop when a < r p . Instead, the radius of the subject drop would be adjusted with equation (25), 
derived in the next section; therefore, equation (20) is of no use for such a situation. 

3.2 Turbulence Characteristic Length and Time Scales 

In the work of Huh et al., the parent drops were assumed to carry homogeneous isotropic turbu- 
lence starting at the injection nozzle exit. 27 Assuming further that no additional internal turbulence is 
generated, the authors were able to obtain an analytical solution of the turbulence scale through the use 
of the well-known k-e turbulence model. Details of this derivation can be found in appendix A. The 
turbulence time scale can be expressed as a function of time as 

r r = r r ° + 0.0828 1. (21) 

The time ( t ) is counted from the time that the parent drop leaves the injection nozzle exit. The initial 
turbulence time scale ( r® ) is evaluated from the initial turbulent kinetic energy ( Ay 0 ) and its correspond- 
ing dissipation rate ( e ® ) at the injector exit: 


, 0 

r° - C — 
1 ~ m o 

£ t 


( 22 ) 


where 
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The velocity (17) is a liquid velocity at the injection nozzle that has the length (L) and the diameter (D). 
The discharge coefficient, the loss coefficient due to the sharpness of the nozzle entrance, and the 
downstream-to-upstream contraction area ratio of the injection nozzle are represented by Q/, K c , and s, 
respectively. Formulations of k r ° and ef are derived in appendix A. 

The turbulence characteristic length scale ( L t ) can be derived in the same way as r t . Huh et al. 
formulated L t as a function of the time and the initial turbulence scales by using the analytical solution 
of the k-e turbulence model: 


L t = L) 


0 (, 0.0828? 


0.457 


1 + - 


(23) 


Similarly to evaluating rf , the initial turbulence length scale ( ) is obtained from the turbulent 

kinetic energy and its dissipation rate as 




k ?r- 


(24) 


Again, details of derivations for the turbulence length scale can be found in appendix A. 

3.3 Product Drop Size and Velocity 

Results from studies of turbulent liquid jets, performed by Tseng et al., show that drops produced 
from the onset breakup regime have sizes in the same order of magnitude of the turbulence length scale. 
On the other hand, with consideration of the surface wave perturbation, Reitz set the product drop 
radius, shown in equation (9), to be proportional to the wavelength associated with the fastest growth 
rate. The above relations suggest that it is reasonable to formulate the radius of the product drops with 
the characteristic length scales for both the wave perturbation and turbulence phenomenon. 

For the present model, the reciprocal of the product drop radius is expressed by the sum of the 
reciprocals of the length scales associated with surface wave instability and turbulence motion, with the 
inclusion of the respective weighting factors as follows: 


-w w_ . 


r P = 


r w r t 




r t + c t r w 


(25) 


The radius of the product drop is r p . The radius (r w ) associated with the wave motion can be determined 
from the Reitz model, as shown in equation (9). The weighting coefficients (c w ) and (c t ) are determined 
by the kinetic energy ratio of the turbulence motion and wave perturbation. These are defined later in 
this section. 
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The value of r t is estimated from the formulations proposed by Huh et al. The authors hypothe- 
sized that the drop size could be represented by a probability density function (PDF) proportional to the 
ratio of the turbulence energy spectrum and the atomization time scale. The notion for this representa- 
tion was that an eddy motion with larger turbulent kinetic energy and a shorter atomization time most 
likely caused the drop breakup to occur more frequently than those containing the lower energy level 
and the longer atomization time. By assigning the wave number of the turbulence energy spectrum as 
the inverse of the product drop diameter, it can be shown that 


f 

j cl 



r 

J a 


x[k e x) 
1 + {k e x) 

{ k e x ) 

1 + {k e x)~ 


-2 1 11 / 6 


lll/6 


dx 


dx 


(26) 


where 

k e = 0.75/L, 

d max = Maximum possible diameter of the product drop 
d min = Smallest possible diameter of the product drop. 

The detailed derivation of the equation above can be found in appendix A. The weighting coeffi- 
cient (c t ) used in equation (25) is simply formulated with the kinetic energy terms related to the surface 
wave and flow turbulence levels as 


c t = 


-‘kt 


Ekt + E 


kw 


(27) 


where 




with L w = A and r w = a/ AQ . 


The relation of c t and the kinetic energy ratio E^ / Ej iW is graphically displayed in figure 4. The 
expression of c t , as shown in equation (27), is obeyed to three constraints: (1) When E^ t «: £/ av , c t 
approaches zero; (2) when E kt » £/ cu ,, c t becomes 1; and (3) if E kt = E kw , c t is set equal to 0.5. These 
constraints impose a balanced contribution of the wave perturbation and turbulence effect on the liquid 
jet breakup and also provide a smooth transition from a given physical dominance to another one. 
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Figure 4. Variation of c t in terms of x w !r t . 


Since the weighting coefficients represent degrees of the contributions of a particular physical 
phenomenon to the overall one, it is logical to assign a unity relation between them as 


c w Q 1 • 


(28) 


With this relation, it is easy to recognize from equation (25) that for the constraint of E^ t » Ef . w , r p 
would become equal to r t , and vice versa, r p would approach r w if E / cw . In other words, equa- 

tions (25) and (27) ensure that the condition with a larger kinetic energy would have a more significant 
effect on the atomization process. 


Since the framework of the present model closely follows that of the Reitz model, the velocity 
and properties of the product drops in the present model are formulated in the same manner as the one 
shown in his model. In each breakup event, product drops are given the same fluid properties and 
physical location as the parent ones. For the velocity, the product drops would carry the same veloc- 
ity (V) of the parent along with two additional normal velocity components (v) and (w) in reference 
to the trajectory of the parent drop. They are represented by 


v = |V|tan(0/2)sin</> 
w = |V|tan(0/2)cos</> , 


(29) 


where 


tan(0/2) = A X AQ/U. 

A random parameter (cp) is chosen on the interval (0, 2ji). According to Reitz and Bracco, the constant 
(Ai) depends on the nozzle design. 40 Its value is set equal to A\ = 0.188 for the present study. 
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3.4 Turbulent Kinetic Energy of Product Drop 

Along with the inclusion of turbulence effect on the primary atomization process, the subject 
phenomenon is also considered in the secondary droplet breakup model, which is discussed in section 4. 
At any rate, this new model requires the initial velocity fluctuation quantity of the product drops right 
after their formation. This quantity can be obtained by examining the energy conservation during the 
primary breakup process. 

When denoting cibe and a a f as the radii of the parent drop before and after its liquid stripping, 
respectively, the amount of the mass stripped from the parent drop is formulated as 


3 3 

fflbr = ~ npi | U[ } e ~ a af 


(30) 


During this breakup process, the change in the individual energy forms of the parent drop is 
estimated as follows: 


Energy due to surface wave motion: ( AE W ( AQ) Z 


Surface tension energy: ( AE sur j ) = 2jrcr|r^ ~ r af 


Turbulent kinetic energy: (AE) Mr ) . = ( — 

\ T t ) 


i \ 

Kinetic energy of motion: ( AE k) par = m br “ 


For the product drop, the energy associated with the surface distortion is negligible, as compared to the 
other energy forms, since the initial disturbance is small. Hence, the energy of the product drop is 
composed of surface tension, turbulence, and kinetic motion as 


N prod 2-Jttp 


Surface tension energy: ( E sur f j 
Kinetic energy of motion: ( E k ) pmd = 


V, 


prod 


where 

Vprod = velocity of product drops 
Nprod = number of product drops 
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The number of product drops ( N proe j ) is obtained from equation (29) while V pr od is calculated from the 
stripped mass (/»/„-) as 

Vprod =Vv“+w 2 , (31) 

Nprod ~ -^3 ■ (32) 

4 TP/ 7 /? 

By equating the change in the energies of the parent drop with the energy of the product drops, the 
turbulent kinetic energy for the product drop is written as 

i^ E w) par + [^ E su>f) + {AEtur) par + {AEk) par ~{( E surf) prod + ( E k) prod ) 

( E 'ur) prod 7) 1 (33) 

n prod 

Equation (33) is used to determine the initial turbulent kinetic energy for the secondary breakup process. 
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4. MODELING OF TURBULENCE EFFECT ON SECONDARY DROPLET BREAKUP 


This section focuses on the formulation of a term to account for the turbulence motion in the 
secondary droplet breakup (TTAB), and the method of estimating the initial turbulence quantities used 
in the TTAB secondary atomization model. The well-known k-e model is employed in the derivation, 
and its initial values are estimated based on the parent drop formation. The turbulence intensity within 
the drops decays as they break up into smaller drops and travel downstream. A criterion of eliminating 
the turbulence consideration in the secondary breakup process is also formulated in this section. 

4.1 Turbulence Term in Secondary Droplet Breakup Model 

Studies of the turbulence effect on the liquid breakup suggest that the turbulence motion tends to 
weaken the surface tension force. 16 In fact, it is well recognized that this surface tension keeps the liquid 
drop from being torn off, while the turbulence within the drop promotes the droplet disintegration proc- 
ess. To account for this behavior, a term (F t ) is introduced to the right-hand side of equation (12) repre- 
senting a force associated with the effect of the turbulence on the droplet breakup as 

= F + F t -kt,- dt, . (34) 

It is assumed that a portion of the internal liquid turbulence energy forms the turbulence 
force (F t ). Another portion of subject energy decays through a normal dissipation process. This proposi- 
tion is based on the experience that turbulence motion behaves like a force participating in the droplet 
surface deformation and accelerating its surface displacement. Hence, it is proposed that the product of 
the turbulence force and the deformation rate of the droplet surface displacement must be related to the 
dissipation rate of the turbulence energy during this breakup process. This relationship is formulated as 

F t t = C t me . (35) 

The dissipation rate of the turbulent kinetic energy per unit mass of individual drops is e and an empiri- 
cal constant representing the proportionality of subject relationship is C t . 

It should be noted that without additional turbulence generation, the internal turbulent kinetic 
energy would decay with time as the drops travel downstream. With an assumption of the turbulence 
within the drop being homogenous and isotropic, the well-known k-e model that describes the relation 
between the turbulent kinetic energy (k) and its dissipation rate (e) can be simplified as follows: 41 

dk , de e 2 . , , __ 

— = -e and — = -Co — , with c P =1.92 . 
dt dt k 

Knowing the initial turbulent kinetic energy and the initial dissipation rate, denoted as ko and eo respec- 
tively, k can be expressed with a function of time as follows: 
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k = 
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(36) 


The derivation of equation (36) is shown in appendix A. It should be noted that ko and eq are determined 
immediately after each droplet breakup process. The expression of £ can now be determined by taking a 
derivative of equation (36) with respect to time as 
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(37) 


The turbulence force, shown in equation (35), can be rearranged and written as 


Fl 

m 


= C t £ o 



(38) 


By substituting equation (38) and the formulations of coefficients defined in equation (12) into equa- 
tion (34), this equation becomes 


£ = C/ A— + C f£o A 
Pl ')> k 0 






Pl r p 
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Pl r p 


(39) 


Again, k is a function of time and is determined from equation (36), whereas C t is based on available 
experimental data. The discussion of estimating C t is presented in section 5. When the distortion dis- 
placement (£) is nondimensionalized by y = , equation (39) can be written in a similar form 

to equation (13) as 
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(40) 


Except for the second term on the right-hand side, equation (40) is the same as the governing equation 
of the TAB model. As stated in section 2, the TAB model has an exact analytical solution. In contrast, 
the exact solution of the present model cannot be obtained since it has a form of a nonlinear and 
nonhomogenous differential equation. Therefore, a second-order central finite difference scheme will be 
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used to solve this equation. The implementation of this model in the numerical calculation is described 
in section 5. 

A similar method, as used in the TAB model, is employed to estimate the postbreakup drop 
size. The SMR (> 32 ) of the new drops is determined from the energy balance between the parent 
and product drops. Derivation of the SMR was presented in section 2. Its formulation is shown here 
for completeness: 


SMR = 


4jlr p O 
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Jt 

~6 


pl 


r " Piy 
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The total energy of the original drop prior to the breakup (E par ) is determined from equation (17). 

4.2 Estimation of Initial Turbulence Quantities 


The secondary breakup process is considered for all drops created from either the primary 
breakup or previous secondary droplet breakup and needs initial values of the turbulent kinetic energy 
and its dissipation rate. These initial turbulence quantities of drops are estimated differently depending 
on the original formation of the drop. 


When a drop is created from the secondary atomization, it is assumed that the turbulence quanti- 
ties of the parent drop are preserved and distributed evenly to its children drops. The turbulent kinetic 
energy (kfa) and corresponding dissipation rate (stn) of the parent drop at the breakup time are determined 
from equations (36) and (37). The time and initial turbulence values used in these equations are refer- 
enced to the time when the considered drop is formed. Formulations employed for this calculation are 


ht = 


M Ce 


(1 




I'/f'-C,) 


(41) 


and 


£ bt = ~ 


£ 0 

,C c 


£ 0 


( k 0 ) 


c, 


-( 1 ~ c e) t bt + { k of l 


n c g /(i-c £ ) 


(42) 


Again, tbt, ko, and sq in equations (41) and (42) are the breakup time and the initial turbulence quantities 
of the parent drop. Based on the conservation of mass, the initial turbulent kinetic energy and dissipation 
rate of the new drop are simply formulated as 
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(43) 
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and 


£ 0 = r ^—^r £ bt • (44) 

SMR 3 

It should be emphasized that ko and sq are denoted for the initial turbulence quantities of a correspond- 
ing drop type. The parameters in equations (41) and (42) belong to the parent drops, while they are for 
the product drops when used in equations (43) and (44). 

When a drop of interest is directly produced from the primary atomization, its initial turbulence 
energy is determined from the energy conservation of the blob drop and its product drops at the breakup 
time. Formulation of the initial turbulent kinetic energy for the product drop, denoted as ( E tur ) pro( j , has 
been derived and already presented in equation (33). Hence, the initial turbulent kinetic energy (ko) per 
unit mass for the drop is represented by 


t »-7 i rfeU’ (45) 

4 JCTpPl 

In their studies of the turbulence droplet breakup, the authors have observed that the turbulence 
characteristic length scale (L f ) has the same order of magnitude as the product drop size. 41 Therefore, the 
product drop diameter is selected from the length scale of the initial turbulence in the present study. The 
initial dissipation rate of the turbulence energy ( £q) is then determined from the relationship of the turbu- 
lence parameters as follows: 
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When L t = 2 r p , so can be written as 
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(46) 


It should be emphasized that equations (45) and (46) are used for estimating the turbulence values of the 
new drops for the primary breakup only. 
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4.3 Criterion of Terminating Turbulence Consideration 


In a normal atomization process, without consideration of droplet coalescence, the liquid drop 
breaks up into smaller drops. As a drop travels downstream its size reduces to the smallest possible 
dimension, and the turbulence activity within the drop also decays. Subsequently leading to conditions 
in which the turbulence approaches its energy cascade regime. Here, the molecular viscosity is effective 
in dissipating the turbulent kinetic energy, and the turbulence scales also reduce to the order of the Kol- 
mogorov scale magnitude. 42 For these circumstances, the effect of turbulence on the droplet atomization 
is no longer necessary since the product drop diameter is either equal to or less than the integral length 
scale, that is approximately equal to four times the turbulence length scale. 41 Obviously, this length scale 
is equal to the Kolmogorov length scale in this limit. Hence, a criterion for eliminating the turbulence 
term in equation (40) is postulated as follows: 
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The bracketed term with its exponent represents the Kolmogorov length scale. It is based on the liquid 
properties and the initial dissipation rate of the turbulence energy of the new drop immediately after 
droplet breakup. When the expression above is valid, the turbulence activity of the new drop is no longer 
effective, and the present model would become the same as the classical TAB model. As seen from the 
results of all test cases in this study, the new drop size reaches the Kolmogorov scale in about two or 
three cycles of breakup steps. 

The proposed Tblob and TTAB models, along with their KH and TAB counterparts, are used 
to simulate various flow test conditions. Their predictions are assessed by comparing the analytical 
results with the measured data. Numerical implementation of these models is described in section 5. 
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5. NUMERICAL IMPLEMENTATION OF ATOMIZATION MODELS 


The first part of this section describes two computer program modules written to evaluate the 
present models. The proposed finite difference schemes of the Tblob and TTAB governing equations are 
formulated here, and modeling simplifications used in the program are stated. Other KH and TAB 
atomization models are used to compare result variations. The last portion of this section covers the 
implementation of these two physical atomization models into an existing CFD program. The discussion 
includes numerical techniques of solving the models, simulated physical domain, and computational grid 
mesh. Sensitivity studies of solution dependency on the computational grid size and selected incremental 
time steps are also presented. 

5.1 Computation of Tblob and TTAB Models for Simple Flow Conditions 

Prior to integration of the proposed models into the existing CFD code, two computer programs 
were written in FORTRAN language to examine the Tblob and TTAB models for simple flow condi- 
tions, in which the liquid jet is injected into a quiescent gaseous ambient environment. Listings of these 
FORTRAN programs are presented in reference 61. The codes track the breakup of an isolated liquid 
droplet as it travels downstream. Fluid properties, initial velocity, and size of the liquid droplet along 
with gaseous conditions such as viscosity, pressure, and temperature are specified and assumed constant 
with time and space. No turbulence in the gas phase is considered for this computation. On the other 
hand, turbulence inside the liquid drop is considered within the present models. Although the program- 
ming is simple and convenient, this framework offers a powerful way to fully assess subject models and 
compare their predictions with the results of the existing atomization models. 

For the primary atomization, the relative velocity between the liquid phase and the local gas flow 
is approximated using the liquid injection velocity. This assumption is reasonable for analyzing this 
breakup regime since this phenomenon occurs mostly near the injector exit. A single blob (parent) drop 
is computed from the injection nozzle exit until it is completely broken up into the children (product) 
drops. The parent drop size changes with time according to equation (20). This atomization process will 
produce the children drops and their size is determined from equation (26). Other formulations for the 
model closure were presented in section 3. 

For the secondary droplet breakup, the TTAB model is used to analyze the atomization process 
of a single drop in a quiescent flow. The governing equation of the TTAB model is solved using a 
second-order finite difference method. Its formulation is shown in the next section. Drop collision and 
coalescence are not included in this calculation; however, they are fully incorporated in the CFD simula- 
tion. In contrast to the primary breakup treatment, a temporal variation of relative velocity between the 
drop and gas phase is computed for the secondary atomization. Hence, the individual force terms in the 
governing equation can be calculated in a more accurate fashion. 

The code for the secondary breakup is composed of two parts. The first part determines the coef- 
ficient (C r ) shown in equation (35). The drop-size correlation formulated from shock tube experiments. 
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conducted by Chou et al., is employed to estimate C t . 44 In the experiments, the breakup of a single liquid 
drop in a flow behind a moving shock is investigated. The average product drop SMR of the entire 
secondary breakup process and the parent drop axial velocity (u p ) are collected and reported in the 
following correlations: 


SMR 2 6 t bt 

r P f * 

and 

u p = 3.754 — • (49) 

t \Pg 

The characteristic breakup time (t ) is defined as 2 r P yjpi/p g /Wg . The time (/) is referenced to the time 
at which the drop is introduced to the flow field, while t Pt is the time when the secondary breakup is 

complete. Results from the subject experiment suggest that t^/t = 5.5 . Therefore, this value is used in 
equation (48) for estimating the SMR of the product drops. The initial radius and relative velocity of the 
parent drop with respect to the gas flow are r p and Wo, respectively. The relative velocity of the parent 
drop was computed from equation (49). The result of C t , presented in section 6, has a value of 0.19. 
Subsequently, this value is selected for the turbulence constant used in this study. 

The second part of this code is used to conduct parametric studies for the test cases listed in 
table 1. Since drops are injected into the ambient gas, the relative velocity is then calculated from the 
following relationship: 
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The drag coefficient ( Cd ) is determined from a model offered by Liu et al. 45 In this model, the drag of a 
deformed droplet is composed of two parts. The first part is the drag of a rigid sphere, while the second 
part involves drag associated with drop surface deformation. The relationship is represented as 
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The Reynolds number (Re) is defined as 2p g Wr p j Droplet distortion displacement (v) is calculated 

from equation (40). When the value of Cd is computed from a previous time step, equation (50) 
becomes an ordinary differential equation and can be conveniently solved. 
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5.2 Implementation of Atomization Models to a Computational Fluid Dynamics Code 


For further evaluation, the present models are incorporated into the CFD-ACE+ program, a 
commercially available software developed by CFD Research Corporation, now owned and marketed by 
the ESI Group. This code, together with its liquid spray module, solves the equations of steady state as 
well as transient flow with chemically reactive fluid dynamics and droplet spray. The governing equa- 
tions and numerical solution method are discussed in the associated software manual and are only 
described briefly here. 50 

The Navier-Stokes equations and renormalized group (RNG) k-e turbulence model for a gas 
phase are solved with a finite volume method. A time marching with a high-order accurate explicit 
scheme is employed in the transient calculation. With an assumption of dilute spray, parcels of drops are 
tracked in a Lagrangian framework. The parcel represents an ensemble of identical drops occupied in the 
same space at a given time. Subsequently, the spray properties at each point are described by statistically 
sampling the spray parcels. The drop parcels exchange mass, momentum, and energy with the gas 
through source terms in the Navier-Stokes equations. 

It should be noted that CFD-ACE+ uses the classical KH and TAB models to predict the primary 
and secondary breakup, respectively. For this study, the existing code is modified to incorporate the 
Tblob and TTAB models. 

5.2.1 Tblob Model in Numerical Simulation 

As stated in section 2, the liquid is assumedly injected in a form of blob parcels containing 
spherical drops with their initial size equal to the nozzle exit dimension. The computational approach for 
the primary atomization proposed by Reitz is adopted in the present study. 22 Extra relations representing 
the turbulence are described here. Initial turbulence quantities are estimated based on equations (21) and 
(22). Values of the turbulence characteristic scales are then calculated from these turbulence values. 
Individual parcels, as described in reference 22, are tracked at each time step, and its drop size is deter- 
mined using equation (20). Next, the mass stripped from the parent drops is computed and is formed into 
a new product drop parcel. Primary breakup of the blob drops continues until their size reduces to the 
dimension of the product drops. 

It is not necessary to create a new parcel in each computational time step. To save computer 
memory and computational time, the stripped mass is accumulated until certain criteria are met. A new 
parcel is then added to the computations and its drop size is calculated from equation (25). The criterion 
of forming a new parcel varies among authors. For instance, Reitz defined the creation of a new parcel 
as when the accumulated mass had exceeded 3 percent of the original parcel mass; while Huh et al. used 
the value of 10 percent. It was found for the test cases of interest, at least, that the computational results 
from theses two criteria are not much different, although the number of parcels in the computational 
domain varies considerably. In CFD-ACE+, a new parcel is created when the number of product drops is 
greater than 20 percent of the ones in the parent. In the present study, the initial turbulent kinetic energy 
of the product drops is estimated from the energy conservation formulated in equation (33) and subse- 
quently equation (45). When this estimated energy is less than 0.1 percent of the kinetic energy of the 
droplet motion, it is assumed that the turbulence motion within the droplet becomes insignificant. No 
turbulence is considered in the secondary breakup for the subject droplet and the classical TAB model is 
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employed. Otherwise, the drops would be considered for secondary breakup with turbulence effects that 
are described in the next part of this section. 

5.2.2 TTAB Model in Numerical Simulation 

While the parcels produced from the primary atomization are continuously tracked in the 
Lagrangian coordinate system, they are also considered for the secondary breakup that is governed by 
equation (40). This equation is solved in the same way as the classical TAB model but it uses the finite 
difference method since its exact analytical solution does not exist. 2 '’ 

To develop the finite difference formulation for equation (40), it is convenient to define a new 
set of coefficients as follows: 
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By substituting the above coefficients, equation (40) becomes 
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In contrast to the TAB model, an exact solution cannot be obtained for equation (40) since it contains 


the extra term C 2 y 1 • Instead, the values of y and y can be approximated with a second-order finite 
central difference scheme in a time step (At) 


y - — 2y" + 

At 2 ' 


-C 3 y n -C 4 


y n+l + y n ~ l 
2 At 


(53) 


29 



To calculate y for the next time step, equation (53) can be rearranged as 

n + 1 _ 2At~C\ + 2At~C2 / n W + 4 — 2At~Cj n + AtC^ - 2 n-1 ,<^ 4 ) 

2 + dtC 4 2 + AtC 4 \‘ > 2 + AtC 4 ' 2 + AtC 4 ' 

For the first time step, at which n = 0, values of y° and y° must be given, and y _1 (that is, the value of y 
at n = -1) must be defined. Furthermore, the first-order backward differencing formulation can be used 
to estimate y ~ l as 


0 -1 

^0 _ y — ^_y _ — 

At 


-1 0 . . - 0 

V =y -Aty . 


Since an initial displacement of droplet deformation and its time derivative are not known for most 
simulation cases, they can be assumed to be equal to zero. Deformation displacement (y) is computed at 
every time step. The secondary breakup assumedly occurs when y approaches 1. The product drop size 
is then calculated from equation (19). At the same time its initial turbulence quantities are estimated 
from equations (43) and (44). 

Again, the new product drops are continuously tracked for possible future secondary atomization 
until they leave the computational domain. When the product drop size is smaller than the Kolmogorov 
length scale, equation (47), the turbulence activity of the new drop is no longer considered. In this case, 
the present model would become the same as the classical TAB model. 

5.2.3 Implementation of Collision/Coalescence 

In the study of the liquid fuel jet breakup, Hiroyasu and Kadota observed that large drops appear 
on the central part of the spray when injecting the liquid with the small spray angle, as defined in 
section 6. 46 Their conclusion was that the large drop results from the collision and coalescence 
processes inside this region. Hence, it is necessary to include these behaviors in the modeling of the 
liquid spray. For this study, a program routine was written to compute the drop collision and coalescence 
model suggested by O’Rourke . 51 Details of the subject model can be found in reference 51. 

A brief review of this model is presented here. 


The collision frequency ( v n ) of droplets contained in any parcel pairs, denoted as parcels 1 and 
2, within each computational cell is examined. This frequency is expressed as 


Vl2 




( 55 ) 


Subscripts 1 and 2 denote drops with a large radius (ri) and a small radius ( 7 - 2 ). The parameter (A3) is 
the number of drops in parcel 2, v denotes the drop velocity vector, and Vol represents the volume of 
the cell. The number of probable collisions (ri) between the droplet pair is calculated from the Poisson 
distribution ( P(n )) as follows: 
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(56) 
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The incremental time (At) is the computational time step. The Poisson distribution ( P(n )) can be stochas- 
tically defined from the uniform distribution in the interval (0,1). When any two drops collide, coales- 
cence between them may or may not occur. Therefore, a collision impact parameter ( b ) is introduced 
to determine these two possibilities. If b is less than a critical value (b cr ), the colliding drops undergo 
velocity changes, while their sizes are unchanged. This phenomenon is called a grazing collision 
reflected by 


b 2 = q(r\ +r 2 ) 2 


b 2 . = (q + r 2 ) 2 min 
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A random number ( q ) is generated on the interval (0,1). When the value of the impact parameter exceeds 
its critical value the drops then coalesce. In this case, n would be the number of drops to be taken away 
from parcel 2. The size and velocity of the new drops must be recalculated according to relations 
described in reference 51. The computation of the collision and coalescence model is carried out 
for all parcels of interest at each time step. 

It should be noted that a considerable coalescence appears when the difference between r\ and r 2 
is large. Furthermore, the collision and coalescence processes are considered for all drops in this study, 
except for the blob drops. It is not logical to have the blob drops included in these processes since they 
are representative of the liquid jet. 

5.2.4 Computational Domain and Grid Mesh 

Eight computational test cases that are related to experiments conducted by various authors have 
been selected for the present work. These test cases will be referred to throughout this TM, and their 
appropriate available measurements and predicted results are presented and compared in section 6. The 
flow conditions of the cases are summarized in table 1. For Tblob, the turbulence scales, initial turbulent 
kinetic energy, and dissipation rate are estimated from equations (22) and (24) by knowing the injection 
nozzle configuration and its flow conditions. Kinetic energy and dissipation rate values are listed in 
table 1. It should be noted that for some test cases the measured data to calculate the initial turbulence 
values are not entirely reported. Subsequently, typical values are estimated and utilized in the calcula- 
tion. Details of determining these initial values are shown in appendix A. 
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Table 1. Test cases and measured data used in the computations. 


Case 

H-1 

H-2 

H-3 

Y-1 

Y-2 

Y-3 

K 

S 

Nozzle diameter (mm) 

0.3 

0.3 

0.3 

0.213 

0.213 

0.213 

0.24 

0.15 

Ambient gas 

Nitrogen 

Nitrogen 

Nitrogen 

Carbon 

Carbon 

Carbon 

Nitrogen 

Nitrogen 





dioxide 

dioxide 

dioxide 



Ambient pressure (MPa) 

1.1 

3 

5 

4.5 

2.5 

0.5 

2.17 

1.5 

Ambient temperature (K) 

298 

298 

298 

298 

298 

298 

298 

289 

Density (kg/m 3 ) 

12.36 

33.7 

56.17 

72.61 

40.34 

8.07 

24.51 

16.84 

Liquid fuel 

Diesel fuel 

Diesel fuel 

Diesel fuel 

Diesel fuel 

Diesel fuel 

Diesel fuel 

Diesel fuel 

Diesel fuel 

Density (kg/m 3 ) 

840 

840 

840 

840 

840 

840 

840 

840 

Viscosity (kg/m s) 

2.9x1 0“ 3 

2.9x1 0" 3 

2.9x1 0- 3 

5x1 0- 3 

5x10- 3 

5x10- 3 

5x10- 3 

2.9x10- 3 

Surface tension (N/m) 

2.05x10-2 

2.05x10- 2 

2.05x10- 2 

2.06x10- 2 

2.06x10-2 

2.06x10-2 

2.06x10-2 

2.05x10-2 

Injection velocity (m/s) 

102 

90.3 

86.41 

185.42 

185.42 

185.42 

133.81 

183 

Initial Turbulence Quantity 

Kinetic energy (m 2 /s 2 ) 

2.88x10 2 

2.26x1 0 2 

3.94X10 1 

1.58x10 3 

1.58x10 3 

1 .58x1 0 3 

5.22x10 2 

4.64x1 0 2 

Dissipation rate of kinetic 

1 .06x1 0 8 

7.34x1 0 7 

1 .22x1 0 7 

1.49x10 9 

1.49x10 9 

1 .49x1 0 9 

3.14x1 0 8 

6.11x10 s 

energy (m 2 /s 3 ) 









Authors 

Hiroyasu and 

Hiroyasu and 

Hiroyasu and 

Yule et al. 

Yule et al. 

Yule et al. 

Koo 

Schneider 


Kodota 

Kodota 

Kodota 






Reference 

46 

46 

46 

47 

47 

47 

48 

49 


Two-dimensional axisymmetric CFD simulations are performed using CFD-ACE+ to compute 
the various cases listed in table 1. A physical domain of 100 x 30 mm in the axial and radial directions, 
respectively, is discretized into a grid mesh of 41 x 25 points in the respective directions. As depicted in 
figure 5, the grid points in this baseline mesh system are clustered near the injector nozzle exit and the 
centerline. The smallest cell of 0.7 x 1 mm in the radial and axial directions contains the injection exit 
port. Boundary conditions used in this numerical simulation are labeled in figure 5. 
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Figure 5. Baseline two-dimensional axisymmetric computational grid mesh. 


The liquid jet, in a form of blob drop parcels, is injected into a quiescent gas at the location of 
x = 0 and r- 0. These drops have the same radius as the nozzle exit. They are injected randomly within a 
prespecified spray angle of the liquid jet at an axial velocity equal to the liquid jet velocity. Although the 
spray angle sensitivity study is not presented here, the results reveal that the variation in the spray angle 
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from 8-18° has only some minor effects on the spray tip penetration. Tip penetration is defined as the 
distance of the parcel that travels the farthest downstream at a given time. With suggestions of various 
authors, the spray angle of 15° is used for this computation. 2427 

It is necessary to point out that only a single parcel each time step is injected at the nozzle exit. 
The number of drops in each initial parcel is determined from the injection mass flow rate, blob drop 
size, and incremental step time. A longer step time would have a larger number of drops per initial 
parcel. Subsequently, a computation with a larger step time would carry fewer parcels as compared to 
one with a smaller step time. However, a large step time may compromise the computational accuracy. 
For this study, a step time (At) of 2.5 x 1CF 6 s is used. It should be noted that At =2.5 x 1(T 6 s is selected 
to have the breakup time of drops at least one order of magnitude larger than the computational step 
time. Small differences are noticed in the computational results for cases of At = 1 x 1 0 6 , 2.5 x 10 6 , 
and 1 xl(T 5 s. 

To investigate the solution dependency on the grid mesh, test case H-l is simulated with four 
grid mesh sizes of 31 x 21, 41 x 25, 61 x 35, and 81x41. Predicted values of the tip penetrations along 
with the measured data are plotted in figure 6. 



Figure 6. Tip penetration predictions with various grid mesh sizes. 


The results indicate that the predicted tip penetrations for all cases are close. In comparison with 
the experimental data, the numerical simulation slightly underpredicts the penetration length in the first 
millisecond and overpredicts afterwards. It should be pointed out that, to some degree, the solution 
dependency on the grid size should be expected, although the present results do not show this trend. 

In a computational study of liquid jet atomization with the dilute spray assumption, Khosla and Crocker 
observed that some solution sensitivity on grid mesh does exist. 52 In addition, a very fine grid system 
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may produce unrealistic results since the dilute spray assumption may no longer be valid. In summary, it 
is reasonable to select the grid mesh of 41 x 25 for this study. 

Computation results of the present models are discussed further in section 6. Predictions of the 
present models and the two classical atomization models are compared and discussed. The study results 
are also appraised against the measured data. 
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6. RESULTS AND DISCUSSIONS 


Prior to implementing the present models into CFD-ACE+, two separate computer routines were 
written to examine the Tblob and TTAB models for simple flow conditions in which the liquid jet is 
injected into a quiescent gaseous ambient environment. The predictions of the present models and the 
two classical atomization models are compared and discussed. The computational results from these 
codes are presented detail. Predictions of the present models are evaluated against the measured data. 
Appropriate comparisons with the existing atomization models and results of the numerical simulations 
are also shown. 


6.1 Primary Liquid Jet Breakup Model (Tblob) 

As stated in sections 2 and 3, the reduction rate of the parent drop size is governed by 
equations (11) and (20) for the KH and Tblob models, respectively. These equations also describe the 
mass stripped from the parent drop as it travels downstream. The term ( C a L t /r t ) on the right-hand side 
of equation (20) represents the contribution of turbulence effects on the rate of change in the parent drop 
size, while the other term (a - C a L w )/r w represents the effects of the surface wave perturbation. The 
values of these terms for test cases H-l and H-3 (table 1) can be found in figure 7. 



Figure 7. Values of the wave motion and turbulence terms on the reduction rate 
of the blob (parent) drop size. 
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In figure 7, the values of these terms are plotted against the relative lifetime of the parent 
drop. This parameter is nondimensionalized by injection velocity and injection nozzle diameter. Results 
for the considered cases indicate that the surface wave perturbation has a considerable effect on the 
reduction rate of the parent drop size. The value of the wave motion term in test case H-3 is approxi- 
mately two orders of magnitude greater than the turbulence term. However, when the injection velocity 
increases from 86.4 m/s (case H-3) to 102 m/s (case H-l) the gap of these two values becomes smaller. 
This suggests that the value of the turbulence term rises at a faster rate than that of the wave motion term 
when increasing the injection velocity. The reduction rate of the parent drop decreases with the increase 
of the velocity. 

It is evident from figure 7 that the reduction rate in the parent drop size is nearly constant with 
time. This trend is due to the fact that the term associated with the surface wave motion is determined 
from the wavelength and its fast growth rate, which are not a function of time. On the other hand, the 
turbulence term is calculated from equations (21) and (23) that include the time terms. Combining these 
two equations when calculating the turbulence term only offers a slight decrease in its value with time. 

In addition, the end of the parent drop breakup occurs in a relatively short time. The results show the 
rate of change in the parent drop size is almost constant. Therefore, the curves of the parent drop 
diameters are nearly straight downward, as shown in figure 8. In all test cases, the Tblob model predicts 
the completion of the parent drop breakup slightly earlier than the KH model. 

It should be pointed out that the Tblob model only describes the primary breakup. However, in 
the actual measurement, it is difficult to separate other physical phenomena such as secondary droplet 
breakup, etc. Therefore, only the intact core length of the liquid jet can legitimately be used to compare 
the measurement and prediction. An available correlation of the intact length that has been widely 
employed in literature is used to compare with the present prediction. 53 The correlation is shown as 



The intact length is denoted as L; et , and the constant ( Cj et ) has a value of 10. In recalling the implemen- 
tation of the KH and Tblob primary breakup models, the liquid jet is injected in the form of blob parcels. 
When traveling downstream, the blob (parent) drop strips its mass to generate the product drops. This 
process is completed when the parent drop size is approximately equal to the size of its product drops. 
Hence, it is reasonable to define the intact length of the liquid jet as the traveling distance of the parent 
drop during this process. Figure 9 shows the nondimensional intact length Lj et j D predicted by the KH 
and Tblob models, as well as the relationship of equation (59) for all test cases plotted against the square 
root of the density ratio Jpi j p g . 


The predicted results from the KH and Tblob models reasonably agree with the correlation 
curve. In general, the data points from the Tblob model are closer to the correlation line in comparison 
to those from the KH model. Furthermore, the Tblob model predicts a shorter intact length than the KH 
model. This prediction is consistent with the measured data trends observed by other authors. 8 ' 53 
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Figure 8. Rate of change in parent drop size predicted by KH and Tblob: (a) Test cases 
H-l through FI-3, and (b) test cases Y-l through Y-3. 
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Figure 9. Comparison of predicted and correlation intact lengths of the liquid jet. 
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In the present model, both the liquid jet surface wave perturbation and the turbulence motion 
play a part in forming the product drops. Their drop size, as shown in equation (25), is composed of the 
two radial length scales r w and r t . They are weighted by the kinetic energies of the respective phenom- 
ena when used to determine the radius of the product drops. The values of subject parameters for test 
case H-3 are plotted in figure 10. The curve showing the kinetic energy of the surface wave motion 
(fig. 10(a)) is very much constant with time because the surface wavelength and its corresponding 
growth rate are not a function of time. On the other hand, kinetic energy associated with turbulence is 
initially at a high level due to the turbulence developed at the injector nozzle exit. However, its value 
rapidly decreases through the dissipation process and approaches the level of the surface wave motion. 
This trend has led to define the weighting coefficient curve, equation (27), as that which starts at the 
value near unity and then quickly reduces its value to one half at the end of the product drop generation. 

Figure 10(b) presents the radial length scales formulated from the two considered motions. Using 
the same previous argument of the surface wavelength being constant, it is easy to recognize from equa- 
tion (9) that the radial length scale (r w ) associated with this motion is more or less constant with time. 

In contrast, the scale (r t ) referring to equation (26) increases with time, since the turbulence length 
scale ( L t ) shown in equation (23) also increases with time. When applying the weighting coefficient (Q), 
shown in figure 10(b), in the determination of the product drop size, the resultant value of r p initially 
rises and then gradually drops due to the reduction of C t . This trend can be interpreted such that the 
high-turbulence intensity at the initial primary breakup stage controls the formation of the product drop 
size. Turbulence dissipates as the parent drop travels downstream. This drop formation process is then 
gradually dominated by the surface wave perturbation. For test case H-3, as determined from the Reitz 
model, r w has a value less than 1 am; while r t , formulated by Huh et al., is approximately 20 tun. 2 Due 
to neglecting the drop collision and coalescence effects, the value of r t is more representative of the drop 
size found in the measured data. The same value is also predicted by the Tblob model at the initial 


38 



CD 

to' 


cn 

O 

o 

CD 


O 

CD 

3 



tO 

o 

o 

CD 


Figure 10. Parameters used to determine product drop size: (a) Kinetic energies and weighting 
coefficient and (b) radial length scales and weighting coefficient. 


product drop formation. However, the drop size continues to decrease to the level predicted by the 
KH model. As the parent drop strips its mass and consequently reduces its size when traveling down- 
stream without coalescence, it is logical to expect that the product drop size should also decrease 
in this process. 


39 


It should be noted that the same trends of all parameters describing test case H-l are also 
observed in all the test cases shown in table 1. The prediction of product drop sizes for test cases H-l 
to Y-3 are displayed in figure 1 1. 




Figure 11. Product drop size predicted by KH and Tblob models: (a) Cases H-l through H-3 
and (b) cases Y-l through Y-3. 
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The results indicate that when liquid is injected at a higher velocity (test cases H-l and Y-l) the 
product drops are initially formed with a larger size. This product drop formation process is prolonged 
in comparison to cases of lower injection velocity (test cases H-3 and Y-3). Similar to the observation 
from figure 10, the results shown in figure 1 1 reveal that the product drop size predicted from the Tblob 
model is approximately one order of magnitude larger than the one from the KH model. Again, by 
ignoring other effects, the drop size estimated from Tblob is consistent with the experimental data. It 
should be restated that the gaseous flow field as well as drop collision and coalescence are implemented 
in the CFD simulation. An effect of the initial turbulent kinetic energy on the primary atomization proc- 
ess is also investigated in this study. The purpose here is to assess the sensitivity of the primary atomiza- 
tion due to the different inlet turbulence levels. This is achieved by using various levels of the initial 
turbulence kinetic energy (ko) while keeping its dissipation rate (t:o) constant. Figure 12 portrays the 
change in the liquid jet intact length (Lj et ) due to the variation of initial turbulent kinetic energy (ho). 

Percentage values of the change in Lj et and ko, shown in figure 12, are referenced to test case 
H-3. The results indicate that the initial energy has a considerable influence on the jet breakup. The 
intact length becomes shorter when increasing the value of ko- It implies that the higher the initial turbu- 
lent energy the faster the primary atomization process accelerates, which leads to an early completion of 
the breakup. Phinney and Tseng et al. also drew similar conclusions from their experimental observa- 
tions. 8 ' 14 It should be noted that the curve in figure 12 approaches asymptomatically to the ko coordinate 
for the higher value of ko- This trend indicates that the initial turbulent energy has a lesser effect on the 
change in the intact length when this energy further increases. 



Figure 12. Change in the liquid jet intact length due to variation of the initial 
turbulent kinetic energy for test case H-3. 
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As the product drops from the primary liquid jet breakup continue traveling downstream, they 
are exposed to external as well as internal forces. Subsequently, they may undergo additional breakup 
cycles. This phenomenon is driven by the secondary breakup mechanism. 

6.2 Secondary Droplet Breakup Model (TTAB) 

A computer routine was written to assess the secondary droplet breakup model. The code tracks 
the breakup of an isolated liquid droplet traveling in a gaseous medium. The model is computed with a 
temporal variation of relative velocity between the drop and the gas medium. Hence, the force terms in 
the governing equation of the drop breakup process can be predicted in a more realistic manner. The 
relative velocity relationship used in this computation was discussed in section 5. 

6.2.1 Evaluation of Turbulence Force Constant (C f ) 

In equation (35), the coefficient (Q) represents the proportionality between the dissipation rate 
of the turbulent kinetic energy and the turbulence energy participating in droplet deformation. For the 
present study, the value of C t is estimated from measurements of drop breakup provided by Chou et al. 44 
In their shock tube experiment, the breakup of a single liquid drop in a flow behind a moving shock 
wave was investigated. The gas in the driven section flows at a velocity of -80.8 m/s with a density of 
1.48 kg/m 3 . Results from the subject experiment suggest that //„//' is =5.5, therefore, this value is used in 
equation (48) for estimating the SMR of the product drops. For the coding of the secondary breakup 
models (TAB and TTAB), the relative velocity of the parent drop was computed from equation (49). 

Due to the lack of measured droplet turbulence quantities, the initial turbulent kinetic energy (ko) used in 
TTAB is estimated from the parent drop velocity, equation (40). Based on numerical simulations of pri- 
mary atomization for the cases listed in table 1, the ratio of the fluctuation velocity within the drop and 
the drop velocity ranges from 0.09-0.11. Now the ratio of 0.1 is selected for determining the initial 
droplet turbulence energy. Since Chou et al. and Liang et al. indicated that the drop starts the breakup at 
tit =1.5, the initial drop velocity is predicted for this particular time. 44 ' 54 The initial dissipation rate (so), 
required in equation (40), is estimated with the relation shown in equation (46). A single water drop of 
590 jam in diameter was analyzed using density, viscosity, and surface tension constants of 997 kg/m 3 , 
8.94 x 1CF 4 kg/ms, and 70.8 x 10 3 N/m, respectively. The product drop size for this case at various C t 
values is plotted in figure 13. 

Since correlation of the experimental data is not related to C t at all, its value at SMR/r /; = 0. 1 35 
remains constant, as shown in figure 13. Note that the TAB model does not carry the turbulence term so 
results from this model would also not vary with C t . Furthermore, TAB predicts the product drop size at 
SMR/r^ = 0.145, which is roughly 7.5 percent larger than the measured value. The solid line curve in 
figure 1 1 displays the results of the present model with a variation of C t . As would be expected, the 
TTAB model reproduces the same results of the TAB model when C t = 0. With an increase in the value 
of C f , TTAB predicts a smaller product drop size. In other words, results of the present model suggest 
that the secondary atomization process would create a smaller drop size when stronger turbulence exists 
within the parent drop. It is evident from the plot that the predicted SMR of the product drop matches 
the measured data for C t = 0.19. Consequently, this value is selected for the turbulence constant used in 
this study. 
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Cf 

Figure 13. Variation of product drop size due to turbulence force coefficient C t . 


6.2.2 Assessment of TTAB Model 

For this particular study, the relative velocity between the drop and the gaseous environment is 
calculated from equation (50). Since the initial turbulence quantities are required for the TTAB model, it 
is worthwhile to conduct a sensitivity study of the turbulent kinetic energy (ko) on the predicted results. 

Variations of the product drop size and drop breakup time at different kg / 0.5 mW^ values for test case 
H-3 in table 1 are plotted in figure 14. 



Figure 14. Product drop size and breakup time of various initial turbulent kinetic energies. 
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The x coordinate represents the ratio of the initial turbulence energy (Icq) and droplet motion 
kinetic energy l/ 2mW^ . Curves of both the product drop size and the breakup time show that these 
quantities are decreasing almost linearly with an increase of the initial turbulence energy, except for 
ko~0. These results suggest that the level of initial turbulence within the liquid has a strong effect on the 
droplet breakup. It is of interest to note that the secondary atomization process generates smaller drops 
with the inclusion of turbulence. Also, it takes a shorter time to break up the parent drop when turbu- 
lence is considered. 

When implementing the secondary atomization model with the CFD code, initial turbulence 
quantities are estimated from the parent drop formation process, as discussed earlier in section 4. As 

the CFD results have shown, kg / 0.5 mW^ ranges from 0.09 to 0. 1 1 for test case H-3 when the drop is 

formed during the primary atomization. Hence, it is reasonable to assign kg/ 0.5 itiWq =0.1 for the 
current study. 

From the governing equation (40), droplet distortion motion is driven by the interaction of aero- 
dynamic, turbulence, surface tension, and viscosity events. To show the variation of these quantities 
with time, their nondimensional values are plotted with time in figure 15. 



Figure 15. Forces generating droplet distortion based on TTAB. 


This plot depicts that the aerodynamic force has a strong effect on the droplet breakup process, 
which is expected, since the flow conditions for test case H-3 relate to a Reynolds number 

(2 Pgi'pWo / Hg) and Weber number [lp„r p W^ /crj of 3.8 x 10’ and 80.7, respectively. Also, the time 

derivative of the distortion displacement (y) of figure 14 starts at a small value and increases with time. 
Since the turbulence force is inversely proportional to the y derivative, this force consequently is very 
high at the initial time and rapidly changes to a smaller value. In contrast, the surface tension and 
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viscosity forces are initially small and change to large values with time because they are directly propor- 
tional to y andy , respectively. According to the breakup model for the case without any turbulence, the 
drop distortion motion is created by the difference of forces between the aerodynamics and the sum of 
the surface tension and viscosity. Obviously this process is accelerated when turbulence is considered. 

A comparison of y and y computed from TAB and TTAB is presented in figure 16. 

The distortion displacement (y) will increase with time when these forces act on the drop. More- 
over, an imbalance of these forces leads to an increase in the time derivative of this displacement ( y ). 
When y reaches the value of one, the parent drop produces smaller drops. Figure 16 also reflects a con- 
sistency with the results observed in figure 15. Because of an additional turbulence term considered in 
the TTAB model, it predicts higher values of y and y at a given time in comparison to the results of 
TAB. Subsequently, the breakup time determined from the TTAB model is shorter than the TAB model. 

Although no measured data of the secondary breakup are available for the test cases, listed in 
table 1, they are still computed with both the TAB and TTAB models to assess the variation in their 
results. Figure 17 shows a comparison of the predictions from TAB and TTAB. 



Figure 16. Comparison of the droplet displacement (y) and its 
derivative (y ) predicted by TAB and TTAB. 
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Figure 17. Difference in product drop size and breakup time 
predicted by TAB and TTAB. 
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The z coordinate of the plot shows the ratio of aerodynamic and turbulence forces. This parame- 
ter is selected for plotting since the aerodynamic force is generally dominant in the droplet breakup 
process. The right and left vertical coordinates represent the percentage differences in the product drop 
size and breakup time, respectively, as predicted from TAB and the present model. The negative scales 
of these two coordinates indicate that the present model predicts smaller values. The results also suggest 
that the large variation of the two predictions appears at the low aerodynamic/turbulence force ratio, 
where turbulence action becomes significant in comparison to aerodynamic force. However, the varia- 
tion becomes smaller with an increase of this ratio value, where the aerodynamic force is dominant. This 
trend suggests that the large variation between the two predictions results from the large turbulence 
force. When the aerodynamic force is significantly higher it plays a key role in the droplet breakup. For 
this condition, the results from TAB and TTAB are not much different, regardless of how high the value 
is of the turbulence force. 

For further evaluation, the proposed models are implemented into the CFD-ACE+ software and 
numerical simulations of the test conditions are performed. These results are presented in section 6.3. 

6.3 Computational Fluid Dynamics Numerical Simulation 

As described in section 5, the Tblob and TTAB models which represent the primary and secon- 
dary atomization processes, respectively, are implemented in CFD-ACE+ for further evaluation. 
Numerical simulations of several test cases using this code are performed and compared with available 
measured data. Also, variations in the predictions from the two classical models KH and TAB, which 
already exist in this code, and the present models are examined. 
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In this study a two-dimensional axisymmetric flow is simulated with a cylindrical liquid jet 
located at the centerline. It should be noted that the axisymmetric conditions, which are imposed here, 
are for mean values associated with the Reynolds averaged Navier-Stokes (RANS) formulation. 
Although instantaneous turbulence quantities are three-dimensional in nature, the k-e turbulence model 
that is used as the baseline model for this study does not account for nonisotropic effects. However, a 
cylindrical liquid jet dominated by the effective shear stresses and the geometric factor associated with 
the stress/strain relationship in this turbulence model would account for the three-dimensional spreading 
in a statistical sense. Time-transient liquid jet breakup and ambient gas behavior are computed for an 
incremental time step (At) of 2.5 x l(T 6 s. The spray trajectory is tracked in the Lagrangian coordinate 
system. A different time step that is normally smaller than At is applied for this coordinate system. 
CFD-ACE+ recalculates this time step continuously based on preassigned maximum allowable changes 
in drop velocity and diameter. The Lagrangian tracking time step, used in the simulations for this study, 
ranges from 1/3 to 1/10 of At. It should be noted that a solution sensitivity study for various numerical 
parameters, such as grid mesh size, At values, etc., was performed prior to obtaining the results of this 
investigation, as reported in section 5. 

Spray tip penetrations from the exit port into the ambient gas for test cases H-l and H-2 are 
plotted in figure 18 as a function of time. Predictions from both the Tblob/TTAB and KH/TAB are simi- 
lar. Nearly straight lines of the tip penetration curves at the beginning indicate that the tip velocity of the 
liquid jets remains almost constant. Once the jet is fully converted into droplets, the surface area of the 
spray nose is enlarged and its speed is subsequently reduced, due to greater aerodynamic force and drag. 
Hence, the penetration curve has a smaller slope at the later time. Furthermore, the results also 
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reveal that the spray penetrates at a faster rate into a lower pressure gas. In measurement comparison, tip 
penetration in both cases is slightly underpredicted at the beginning of the spray with analytical models; 
however, they are somewhat overpredicted at the later time. Overall, the computational results show 
close agreement with the experimental data. 

Variations in the spray angle due to different ambient gas pressures are depicted in figure 19 for 
test cases H-l through H-3. This angle is a global parameter used to characterize the radial expansion 
of the spray. For an axisymmetric case, as shown in figure 19, the angle is formed by two opposite tan- 
gent lines drawn along the spray envelope. The results portray higher backpressure with a larger radial 
expansion of the spray. Further examination of the atomization modeling suggests that the radial product 
drop velocity, as estimated from equation (29) of the primary breakup process, is related to the wave- 
length and its maximum growth rate. The values of these parameters, however, are dependent on the 
liquid-to-gas density ratio. Consequently, higher gas pressure would result in a larger radial-to-axial 
velocity ratio and a wider spray angle. In comparison with the results of the KH/TAB models, the pre- 
sent models predict a more favorable, larger spray angle. As stated in section 6.1, the Tblob model esti- 
mates considerably larger product drop sizes than the KH model, which leads to the prediction of a 
stronger radial momentum. Consequently, liquid drops are spread outward more in the Tblob model 
when compared to the predictions provided by the KH models. As shown in figure 19, a better agree- 
ment with the measured data is obtained with the present models. 



Figure 19. Variations in predicted and measured spray angles at different backpressures. 


Computed spray shapes and tip penetrations from the present models for three different back- 
pressure conditions after 2.5 m/s of liquid injection are displayed in figure 20. It should be noted that the 
injection velocity is higher for a lower backpressure condition. Hence, the spray tip penetrates faster into 
a lower pressure environment, although the drag increases on the periphery of the spray contour. This 
phenomenon leads to parcels containing small drops being dragged into the wake behind the tip region, 
which can be clearly observed in figure 20(a) and (b). It is also interesting to note that larger drops 
appear on the spray tip as well as on the peripheral regions near the tip. Generally, large drops carry 
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a stronger momentum that decay less than their smaller drop counterparts. Hence, the large drops 
tend to travel downstream faster than the smaller drops. In addition, large drops found in these regions 
may reflect the effect of drop coalescence, which is more pronounced in an elevated gas pressure. As 
figure 20(c) illustrates, a larger drop size can be observed in the tip region at higher ambient pressures. 
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Figure 20. Variations in spray shape, tip penetration, and drop size at t = 2.5 m/s due 
to different backpressures: (a) 1.1 MPa, (b) 3 MPa, and (c) 5 MPa. 


Results collected from numerical simulations of test case K are shown in table 1. This test case 
was selected to evaluate the atomization models because of its available measured drop size data. In the 
experiment, Koo measured diesel sprays with a phase/doppler particle analyzer (PDPA) system. 48 In 
general, a PDPA can only recognize a defined range of droplet diameters due to digital processing 
hardware limitations, as described in the work of Rhys. 55 Normally, the ratio of the largest observable 
droplet to the smallest observable droplet cannot exceed 35. Regarding this experiment, the instrument 
was set to measure drop sizes from 1.8 to 246 um. In this analysis, the same grid mesh and physical 
domain were used as displayed in figure 5. It should be noted that the liquid fuel injection was driven 
by a cyclic pumping system for which the mass injection rate is no longer constant. To more realistically 
simulate the inlet flow conditions, the transient injection velocity profile, reported by Koo and shown in 
figure 21, was utilized in the computations. 

Similar to the previous computational cases, this test case also simulated the transient flow 
field conditions with At = 2.5 x 10 6 s. It should be noted that the present simulation utilizes the RANS 
approach to obtain the mean velocity field of the gas phase. For statistically unstationary flows, which 
were investigated experimentally by Koo in this case, the interpretation of the averaging process is 
ensemble averaging. 48,56 Thus, the RANS formulation is justified for the statistically transient case 
as indicated by the injector inlet condition of figure 21. 
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Figure 21. Fuel injection velocity used in simulation for test case K. 


The predicted spray tip penetration and the measured data are presented in figure 22. Similar to 
the computational results obtained from test cases H-l and H-2, the predictions of this test case, with 
both the KH/TAB and the present models, are slightly lower than the measurement at the beginning of 
the injection; however, the predicted curves approach the experimental data at a later time. It is noted 
that Koo reported the spray penetration length from the beginning of the injection to 1.5 m/s while the 
computation continues to 3 m/s. In general, both the predicted tip penetrations are similar and their val- 
ues are in reasonable agreement with the experimental data. 



Figure 22. Predicted and measured spray tip penetration for test case K. 
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Figure 23 displays the simulated spray shapes along with photographs recorded at several time 
intervals in the experiment. No significant variations appear between the computational results produced 
for the two sets of atomization models, KH/TAB and Tblob/TTAB. In comparison with the photographs, 
the simulated spray penetrates into the ambient gas with a slightly shorter distance than the actual dis- 
tance depicted initially from the photographs; however, the predictions improve significantly afterwards. 
These results also reflect the same assessment of the tip penetration comparison shown previously in 
figure 22. It is also interesting to observe the degree of the spray being spread laterally and its tip shape. 
All the models predict a similar radial expansion of the spray for all considered time intervals and they 
compare reasonably well with their photographic counterparts. However, the numerical simulations dis- 
play the tip profile more like a conical shape with parcels containing large drops rather than a round tip 
as shown in the photographs. Since the momentum of these large drops is degrading much less than the 
small drops, the large drops move downstream faster than the smaller drops. Therefore, the moving 
motion of the large drops creates the appearance of the conical shape at the tip. An analysis of the com- 
putational cases has been conducted to examine the tip shape variation for different initial spray angles. 
The analysis suggests that the spray tip head changes to a round shape and does not penetrate into the 
ambient gas as deeply when increasing the initial spray angle. Tip variation analysis is not addressed in 
detail here since it is not a main focus in this study. 

Further detailed drop size comparisons between the predictions and the measurements are pre- 
sented in figures 24 and 25. Figure 24(a) reveals the drop sizes recorded at three axial locations along 
the centerline of the spray. The results indicate that the jet, with an initial diameter of 240 pm, breaks up 
into relatively large drops along the centerline. Their product drop diameters range from 100 to 200 pm 
at the stations of v= 10 and 20 mm, as shown in figures 24(b) and (c). These results suggest that the pri- 
mary breakup process is dominant in these regions. Further downstream, the drops undergo the secon- 
dary atomization that generates smaller drops. Figure 24(d) exhibits drops of 20-50 pm in diameter at 
the x = 60 mm station. It should also be noted that the fuel injection cycle for this test case is complete at 
1.77 m/s; therefore, the drop sizes shown in figures 24(b) and (c) also reflect this end period at which 
the small drops are registered. The small drops recorded after 1.77 m/s are formed because of the small 
injection rate at the tail end of the injection cycle, as well as the small drops lagging behind in the spray 
motion. Generally, predictions from the present models are in good agreement with the measured data. 

In contrast, the KH and TAB models underpredict the drop size for stations x = 10 and 20 mm; however, 
their predictions are improved significantly at the x=60 mm station. This observation is consistent with 
earlier results presented in section 6.1. The results also show that the KH model generally predicts 
smaller drop sizes than the Tblob model. 


51 





Figure 23. Simulations and photographs of sprays for test case K: (a) At t = 0.203 m/s, 
(b) at t = 0.601 m/s, (c) at t= 1.205 m/s, and (d) at t= 1.8 m/s. 
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Figure 24. Drop-size distribution from the predictions and measurements for test case K: 
(a) Coordinates of locations where drops recorded, (b) drop size at [10,0], 

(c) drop size at [20,0], and (d) drop size at [60,0], 


Figure 25(a) displays the drop sizes at three off centerline locations [x,y] . Relatively large 
drops registered at the [10,0.5] location with diameters ranging from 30 to 180 pm with a random 
order in time as shown in figure 25(b). Smaller drops are observed after the end of the injection cycle at 
t = 1.77 m/s, similar to the results shown in figure 24. The present model predictions capture large drops 
intermittently at certain short time intervals; however, both the KH/TAB and Tblob/TTAB atomization 
models typically project a much smaller drop size than the measurement portrayed. In other words, the 
CFD simulation indicates that the secondary droplet breakup process takes place at this location but not 
the primary breakup. The drop diameter at [10,1.5] and [30,2] varies from 20 to 50 pm, as shown in 
figures 25(c) and 25(d), respectively. This suggests that considerably smaller drops appear on the spray 
periphery, and their presence in this outermost region forms the shape of the spray. Furthermore, since 
the drop sizes at these locations are relatively small, the secondary droplet breakup occurs mostly in 
these regions. This assessment is consistent with the observation of Koo in his measured data analysis. 
He concluded that the small droplets in the periphery of the spray are not those directly injected from the 
nozzle. The small droplets from the nozzle travel much slower than the large droplets that are produced 
in the primary breakup; hence, they cannot maintain the spray shape observed in the experiments 
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Figure 25. Drop-size distribution from the predictions and measurements for test case K: 
(a) Coordinates of locations where drops recorded, (b) drop size at [10,0.5], 
(c) drop size at [10,1.5], and (d) drop size at [30,2]. 


without the appearance of drops formed by the secondary breakup mechanism. It is interesting to 
note that the simulation with the KFI and TAB models has a limited number of parcels and no parcel 
tracked at [10,1.5] and [30,2], as shown in figures 25(c) and 25(d), respectively. Again, this exami- 
nation reconfirms that the spray angle predicted by the classical KFI and TAB models is a little 
smaller than the actual angle, as displayed in figure 18. In the present model, numerical simulation 
captures droplets at [10,1.5] and [30,2] only in a short timespan, but the predicted drop sizes at this 
interval are close to the measured data. 
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7. CONCLUSIONS AND RECOMMENDATIONS 


The main effort of this research is to model the effects of turbulence within a cylindrical liquid 
jet on the atomization process. The present models have been formulated by adding the terms accounting 
for these effects to the widely used KH and TAB atomization models. Predictions from these new mod- 
els have been fully evaluated by comparing the results of the existing models and the available meas- 
urements. This section summarizes the research effort, reports several findings captured in this study, 
and offers recommendations for further research on the subject topic. 

7.1 Summary 

In the course of this study, two existing classical atomization models, the KH instability devel- 
oped by Reitz and the TAB by O’Rourke and Amsden, have been extended to incorporate terms that 
account for the turbulence contribution on the liquid breakup mechanism from a phenomenological 
point of view . 22 ' 23 Subject terms are derived primarily from the physical modeling work of Huh et al. 
and the collective observations of Tseng et al. from their experimental investigations of the turbulent jet 
breakup . 2714 It is important to point out that the present models do preserve the original predictions of 
these two classical models when turbulence is not present. 

For primary atomization, the KH model has been modified to implement the turbulence effect 
using the characteristic scales method associated with individual physical phenomena. In the proposed 
Tblob model, the characteristic length and time scales of surface wave perturbation and those of turbu- 
lence motion are combined, so their contribution to the breakup mechanism is weighted by means of 
kinetic energy. The rationale for this approach is that the larger kinetic energy motion has a stronger 
influence in the liquid breakup process. This has led to modification of the equations predicting the mass 
rate stripped from the blob drop and the product drop size. It has been observed that initial turbulence 
quantities play a key role on jet disintegration. Their values are estimated from the geometry and flow 
conditions of the injection nozzle. 

For the secondary droplet breakup event, an additional force term that represents the turbulence 
effect on droplet surface distortion is incorporated into the TAB governing equations. This force term is 
composed of the deformation rate of the surface distortion, the dissipation rate of the turbulent kinetic 
energy, and an empirical constant. The value of this constant rests on the data collected from the meas- 
urements of drop breakup in a shock tube. Similar to the Tblob primary atomization, the present secon- 
dary droplet breakup model also requires initial turbulence quantities. These values are estimated from 
the energy conservation of the parent drop and product drops from the primary atomization process. 
When a drop is created from the secondary breakup, turbulence quantities are determined simply from 
the turbulence values of their parent drop. 

Two separate computer programs have been written to individually assess these new models 
by comparing the predictions with available experimental data as well as results predicted from their 
classical model counterparts. The proposed atomization models were implemented into an existing 
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CFD-ACE+ code for further evaluation. The drop collision and coalescence model of O’Rourke was 
also incorporated into the CFD-ACE+ code to provide more realistic simulations. The computational 
results were then compared with the measured data and predictions of the original models. Several 
conclusions can be drawn from these investigations and are reported in the following section. 

7.2 Conclusions 

The proposed atomization models have been used to simulate several test cases. The flow condi- 
tions are listed in table 1. Close examination of the results, presented in section 6, has led to the follow- 
ing conclusions: 

• The primary atomization regime predictions generally show that the intact core length of the turbulent 
liquid jet is slightly shorter than without the turbulence consideration. In fact, at least in the test cases 
considered in this study, the combined characteristic scales used in the Tblob model result in a shorter 
jet breakup time, leading to a shorter intact core as compared to the KH model. These predictions are 
consistent with the breakup observations of fully developed turbulent liquid jets by the other 
authors. 8,53 

• When the turbulent liquid jet disintegrates and then forms droplets, as observed by Tseng et al., 
these drops are generally larger than the drops produced by a nonturbulent jet. 14 The computational 
results of the present models also reflect this finding, particularly in test case H-3, where the product 
drop size predicted from the TBlob model is approximately one order of magnitude larger than the 
KH model. 

• In general, turbulence inside the liquid jet affects the primary breakup of the liquid jet. This turbu- 
lence is characterized by fluid properties, gas flow conditions, and initial turbulence quantities ko and 
et). In turn, these initial quantities are dependent on the flow conditions and the geometry of the injec- 
tion nozzle. 

• According to the results from the TTAB model, turbulence inside the parent drop also plays a role in 
the droplet surface wave distortion. Similar to the primary breakup, the level of the turbulence effect 
on the secondary droplet breakup process depends on its initial turbulence values. However, this tur- 
bulence effect is diminished after a few droplet breakup cycles. This is due to turbulence reduction 
from one generation of drops to the next, based purely on the current assumptions of estimating initial 
turbulence quantities. 

• Close examination of the predicted results for several test cases suggest that the values of the forces 
participating in droplet deformation vary significantly at the initial time. Aerodynamic and turbulence 
forces are the strongest among them. However, turbulence force approaches a same order of magni- 
tude as surface tension and viscous damping forces at the breakup time, while the magnitude of the 
aerodynamic force is still maintained. This observation reconfirms that turbulence has a considerable 
influence on the secondary breakup when compared with surface tension and viscosity effects. But, 
aerodynamic force still plays a dominant role in the breakup process also. 
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• In contrast to the characteristics of liquid jet disintegration, the secondary breakup mechanism pro- 
duces small drops with a short breakup time when turbulence is considered. The turbulence force 
term in the present model is constructed to promote droplet surface distortion. That formulation leads 
to predict a smaller drop size and shorter breakup time than the results obtained from the existing 
TAB model. 

• In general, the CFD simulations are able to capture similar spray structures observed in the experi- 
ments. The spray tip penetration lengths as well as the appearance of large drops at the spray tips, as 
predicted from the two sets of atomization models (KH/TAB and Tblob/TTAB), are comparable to 
the measured data. However, the KH and TAB models estimate smaller spray angles for test cases 
H-l through H-3 than the measurements indicate. 

• Drop measurements and predictions from test case K indicate that small drops created from the sec- 
ondary breakup appear in the outermost regions of the spray. This breakup mechanism is responsible 
for the radial expansion of the spray and these drops then form the spray profile. On the other hand, 
drops on the centerline of the spray are relatively large. The appearance of large drops at the spray 
tips may result from drop coalescence. The strong momentum of large drops allows them to travel 
ahead downstream. 

• In summary, for test case K, the CFD simulation using the KH and TAB models tends to underpredict 
the drop size of the spray, while it is able to reasonably capture the liquid spray structure. On the 
other hand, the present models predict comparable drop sizes as seen in the measured data. Overall, 
the predicted results from the present models do agree reasonably with the experimental data. 

7.3 Recommendations 

This research effort has offered enhancements on the widely used KH and TAB atomization 
models by incorporating the turbulence effects. Results from the test cases considered in this investiga- 
tion indicate that the present models make reasonable predictions in comparison with the available 
measured data. Nevertheless, the present form of these models may not always provide satisfactory 
results for all other flow conditions. Hence, to make use of the present models and to further evaluate 
and improve them, the following recommendations are offered: 

• The present models are able to provide global predictions of the atomization process and the spray 
formation for certain flow conditions. Their results are adequate and useful for most practical engi- 
neering applications. Nevertheless, by no means are they able to represent accurate and complete 
descriptions of physical happenings in the liquid atomization process, since they are based on a 
semiempirical approach and a phenomenological point of view. 

• Reasonable results can be obtained from the present models as long as the flow conditions of interest 
are similar to the conditions studied in this research. In addition, the solution is dependent on the grid 
mesh, computational domain size, time step, etc., and should be assessed and anchored, if possible, 
prior to using the computational results. 

• The new models simulating the atomization process are based partially on knowing local flow condi- 
tions. In other words, the predicted results also depend on the local flow solutions that are computed 
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by CFD simulations. Unfortunately, these solutions may vary slightly among CFD numerical 
schemes. Certain care should be taken when implementing these new improved models into a CFD 
code. The predictions should be validated adequately with available measured data. It is recom- 
mended to simulate several test cases provided in table 1. 

• This research provides a basic framework to enhance the atomization models; however, the coeffi- 
cients used in both the original and the present formulations should be reviewed and reevaluated 
when used beyond the investigated flow conditions. Recent results from the investigation by Grover 
et al. 57 may be considered when applying these models to a wide Weber number flow range. Grover 
et al. proposed to improve the atomization models by changing the subject coefficients according 

to the considered flow conditions. 

• The predicted results of the present models are greatly dependent on initial turbulence quantities. 
Unfortunately, the current form of the present models only provides a rough estimate of these initial 
values. To obtain more realistic predictions, a more rigorous estimation of the initial turbulence quan- 
tities for individual injection nozzles of interest should be found. 

• Finally, the present study indicates that, at least for the test cases investigated, the TAB model pre- 
dicts smaller drops than anticipated on the secondary breakup. With the inclusion of the turbulence 
effect, the present TTAB model estimates even smaller drop sizes. The recent study by Lengsfeld et 
al. suggests that the droplet breakup process could be altered by the presence of neighboring droplets. 
However, the present form of the models does not take this effect into account. Further studies on this 
topic are suggested to improve the prediction capability. 58 
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APPENDIX A— FORMULATION OF TURBULENCE QUANTITIES 


This appendix presents detailed derivations of relations and equations involving the turbulence 
consideration in the liquid jet atomization process. These equations are referenced throughout this TM. 
Formulations to determine the initial turbulent kinetic energy and its dissipation rate for use in the 
primary breakup model are shown. Derivations of the turbulence scale relationships and the product 
drop-size scale associated with the turbulence are presented. Finally, estimated turbulence values for the 
test cases used in this study are discussed and listed. 

A.l Initial Turbulent Kinetic Energy 

Consider a liquid jet issued from an injection nozzle that has the length (L) and the diameter (D) 
as shown in figure 26. The outlet-to-inlet contraction area ratio of the nozzle is s. The liquid is injected 
at a velocity of U. 



Figure 26. Configuration and flow conditions of a typical injection nozzle. 


The total pressure drop ( AP tot ) of the injection nozzle is composed of three components: (1) The 
pressure drop associated with the nozzle entrance shape ( APf orm ), (2) the flow acceleration pressure drop 
(AP acc ), and (3) the pressure drop due to the wall shear stress (AP noz ): 

AF \ot = APform + AP acc + AP noz . (60) 

The total pressure drop is expressed in terms of the nozzle discharge coefficient (Q/) as 

( 61 ) 

2 Crf 
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Depending on the shape of the entrance corner, the pressure drop occurs when the fluid enters into the 
nozzle and can be estimated with a loss coefficient ( K c ) as 

APform-Kc^Y- ■ ( 62 ) 

For a contracting nozzle, the fluid is accelerated due to the reduction of the cross-sectional area from the 
nozzle inlet to its outlet. The pressure drop associated with the flow acceleration process is estimated as 

^P acc -(l-s 2 )^y- . (63) 


When substituting the proper terms with equations from (61) to (63) into equation (60), AP noz can be 
calculated as 


AP, 


noz 


= pJC( x . s Ap^1. k pJ!1 


2 C 2 d 


I 2 


(64) 


For a fluid flowing through the injection nozzle, Huh et al. approximated the average turbulent stress by 
setting it equal to the wall shear stress 27 


Pl u 2 xDL = AP noz . (65) 

The parameter (u 2 ) is the square of the turbulent fluctuation velocity. After placing equation (65) into 
equation (64), the initial turbulent kinetic energy ( ) can be calculated as 



W 


&L/D 


~2~ K c 

Q 



( 66 ) 


A.2 Initial Dissipation Rate of Turbulent Kinetic Energy 

To determine the initial dissipation rate of the turbulent kinetic energy, Huh et al. assumed that 
the dissipation rate is proportional to the work done by AP noz on the nozzle volume as follows: 27 


Pi 


xD 2 L 
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: K e AP noz 


nD" 


L 


U_ 

L 


( 67 ) 


The average dissipation rate of the initial turbulent kinetic energy is represented by e z . The 
proportionality constant (K f ) is found equal to 0.23 by comparing it with multidimensional internal 
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nozzle flow calculations. 59 By substituting AP noz in the equation above with equation (64), e® can be 
expressed as 


f £ 2 L 


\-K c -[\-s 2 

cs 


( 68 ) 


A.3 Derivation of Turbulence Scales 


It should be noted that the internal turbulence of the parent drop decays with time as it travels 
downstream. Assuming a homogeneous isotropic turbulence with no internal turbulence generation, the 
well-known k-e turbulence model can be written in terms of the time derivatives of the turbulent kinetic 
energy (k) and its dissipation rate (e) as 


and 



(69) 


de _ £ 2 

— = -C E — 

dt k 


(70) 


The k-s turbulence constant (C f ) is equal to 1.92. The equations above do not account for the effects 
of the parent drop breakup process and collision/coalescence to the drop internal turbulence. Huh et al. 
were able to obtain an analytical solution of the k-e equation model for the individual parent drop as 

follows: When e in equation (70) is replaced with dk/ dt of equation (69), the relation of the derivatives 
between k and e becomes 27 


de 

e 



(71) 


When integrating equation (71) with the initial conditions of k(o) = kf and e(o) = e ® , the solution 
of £is 



t v\ 






Replacing e of equation (71) with equation (72), the derivative of k can be expressed as 




dt . 


(72) 


(73) 
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The turbulent kinetic energy (k) can be obtained by integrating equation (73) over time and with the 
initial condition of k(0) = k E 


k = 


K 


0 \C £ f 1 C e) t + [ k t 




V(l ~C e ) 


The turbulence characteristic time and length scales are defined in terms of k and e as 



L, = C, 




( 74 ) 


(75) 

(76) 


The parameter ( C u ) is another k-e turbulence constant and is set equal to 0.09. r t and L t can be expressed 
as a function of time and initial conditions by performing substitutions of e and k from equations (72) 
and (74) into equations (75) and (76): 


T t~Cf_i(C E l)t + r f 


and 


where 


L r =L? 




1 -c 


£ / 


-C c 


T?-C p %andL?-C f l "' ; 




For C jU = 0.09 and C £ = 1.92, equations (77) and (78) become 

r t = rf + 0.0828t 


(77) 


(78) 


(79) 
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and 


L r =L? 


1 + 


0 . 0828 ? 
Io" 


, 0.457 


( 80 ) 


A.4 Product Drop-Size Scale Associated With Turbulence Motion 

In the present primary breakup model, the product drop size is formulated using two unique 
length scales — surface wave perturbation and turbulence motion. The scale associated with turbulence 
is derived in the remainder of this appendix. 

A.4.1 Governing Equations 

Huh et al., represented the diameter of the product drops associated with turbulence motion by 
a probability density function denoted as P(d t ), where d t is the diameter of the product drop. 27 This 
function is proportional to the turbulence energy spectrum and the reciprocal of the time scale. It is 
based on the notion that an eddy motion with larger kinetic energy and shorter atomization time is likely 
to appear more frequently as an atomized droplet. 


P{d,)-C, 




where d>(d t ) is the turbulence energy spectrum given as 56 




[k/kef 

1 + (k/k e ) 


2 lt/6 


( 81 ) 


( 82 ) 


The wave number (k) is the inverse of the droplet diameter (d t ) and k e is obtained from 


k e 


— with L, 

L e 


u 

0.75 ’ 


where L e is the size of eddies with the maximum energy in the turbulence energy spectrum. The constant 
(' C p ) is determined by the normalization condition 


+ 00 , v 

J P(x)dx = 1 , 


( 83 ) 


where x is a dependent variable of the integral that represents the product drop diameter (dt). 
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A.4.2 Product Drop-Size Scale 

To calculate the diameter of the product drop, the normalization coefficient ( C p ) is evaluated by 
substituting equations (81) and (82) into equation (83) resulting in 


^ p nontax 


(k e x) 


-2 


1 + {kex)~ 


11/6 


dx = 1 . 


(84) 


It should be noted that the integral in equation (84) is not finite due its infinity limits. Therefore, the 
limits of the original integral have been changed from +°° and -oo to d max and d min . For the current 

calculation, d max is set equal to the diameter of the injector nozzle (D) while d m j n is set tolO“ s m. 


Since the product drop size involving the turbulence motion is estimated from the probability 
density function, as shown in equation (83), the associated radius can be written as follows: 


r t = 


df 


~-f ma x p(x)dx = -^-f m 
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x(k e x) 
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dx 


(85) 


Substituting the ratio C p jx t in equation (85) by its equivalent in equation (84), the diameter of the 
product droplets can be calculated from the ratio of the following two integrals: 
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dx 
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dx 


The integrals in the equation above are evaluated with Simpson’s composite algorithm. 60 

A.4.3 Initial Turbulence Quantities 


( 86 ) 


For the Tblob primary breakup model, the initial turbulent kinetic energy and its dissipation rate 
are required. These values, which are estimated from equations (66) and (68), need injection velocity, 
nozzle geometry, and nozzle flow parameters as shown in table 2. This table lists values of subject 
quantities used in the test cases and stated in section 5. The test case numbers are labeled in the same 
way as shown in table 1. Several values in this table are assumed and/or estimated due to unavailable 
measured data in some experiments. These are denoted by a footnote under the table. 
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Table 2. Parameters of test cases for estimating turbulence quantities. 


Case 

H-1 

H-2 

H-3 

Y-1 

Y-2 

Y-3 

K 

S 

Nozzle diameter, D jnj ( mm) 

0.3 

0.3 

0.3 

0.213 

0.213 

0.213 

0.24 

0.15 

Nozzle length, L (mm) 

0.8 a 

0.8 a 

0.8 a 

0.67 

0.67 

0.67 

0.8 

0.8 a 

Discharge coefficient, C d 

0.7 b 

0.7 b 

0.8 b 

0.6 

0.6 

0.6 

0.67 

0.7 b 

Entrance loss coefficient, K c 

0.45 c 

0.45° 

0.45 c 

0.45 c 

0.45 c 

0.45 c 

0.45 c 

0.45° 

Injection throat/exit area ratio 

0 d 

0 d 

0 d 

0.01 

0.01 

0.01 

0.01 

0 d 

Injection velocity (m/s) 

115.8 

102.54 

86.41 

185.42 

185.42 

185.42 

133.81 

183.0021 

Initial Turbulence Quantity 

Kinetic energy (m 2 /s 2 ) 

2.88 xIO 2 

2.26 xIO 2 

3. 94x10' 

1.58x10® 

1.58x10® 

1.58x10® 

5.22 xIO 2 

4.64x10 2 

Dissipation rate of kinetic 
energy (m 2 /s 3 ) 

1.06x10® 

7.34 xIO 7 

1.22x10 7 

1.49x10 9 

1.49x10 9 

1.49x10 s 

3.14x10® 

6.11x10® 

Authors 

Hiroyasu and 
Kadota 

Hiroyasu and 
Kadota 

Hiroyasu and 
Kadota 

Yule et al. 

Yule et al. 

Yule et al. 

Koo 

Schneider 

Reference 

46 

46 

46 

47 

47 

47 

48 

49 


a Assumed value due to the lack of measured data. 
b Estimated based on given fuel conditions at the nozzle exit. 
c Assuming that the injection nozzle has a sharp entrance corner. 
d For the straight injection tube, without throat. 
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